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Fortran  Programs  for  the  Calculation  of  Wigner  3 j , 
6j,  and  9j  Coefficients  for  Angular  Momenta  ^  80 

R.  S.  Caswell  and  L.  C.  Maximon 

Fortran  II  and  Fortran  IV  programs  are  given  for 
the  calculation  of  Wigner  3 j ,  6j,  and  9j  coefficients 
containing  individual  angular  momenta  less  than  or 
equal  to  80.   The  large  numbers,  resulting  from  the 
factorials  which  appear  in  the  expressions  for  these 
coefficients,  are  handled  by  taking  logarithms  and 
performing  most  of  the  pertinent  arithmetic  with  the 
logarithms.   The  alternating  series  involved  in  the 
expressions  for  the  3j  and  6j  coefficients  are  summed 
arithmetically  in  double  precision.   We  present  two 
versions  of  the  basic  program  for  the  calculation  of 
these  coefficients  (each  version  being  given  both  in 
Fortran  II  and  in  Fortran  IV).   All  programs  are  called 
as  single-precision  Fortran  functions.   The  difference 
between  Versions  I  and  II  is  in  the  internal  structure 
of  these  programs.   In  Version  I,  only  the  series  summa- 
tion is  in  double  precision,  resulting  in  a  faster  pro- 
gram which  occupies  less  memory  in  the  computer,  the 
computation  time  being  of  the  order  of  0.001  to  0.01 
sec  per  6j  and  0.1  to  1  sec  per  9j.   Version  II,  being 
entirely  in  double  precision  internally,  is  much  more 
accurate  (the  errors  being  of  the  order  of  one  percent 
of  the  errors  in  Version  I),  but  the  computation  time 
for  a  given  coefficient  is  approximately  50%  greater 
than  in  Version  I.   Information  relative  to  the  accuracy 
and  memory  requirements  of  the  two  versions  of  the  pro- 
gram is  provided,  and  Fortran  II  and  IV  lists  of  the 
programs  and  calling  instructions  are  given. 

Key  Words:   Fortran  program,  Wigner  coefficients, 
quantum  theory  of  angular  momentum, 
3j  coefficient,  6j  coefficient, 
9j  coefficient,  Clebsch-Gordan  coefficients 

1.   Introduction 
Existing  programs  for  calculation  of  Wigner  3 j ,  6j,  and  9j  coef- 
ficients of  which  we  are  aware  either  (1)  do  not  work  at  all  for  large 
values  of  the  arguments,  or  (2)  produce  significant  errors  (for  example, 

~  20  percent)  when  used  for  large  values  of  the  arguments,  or  (3)  exist 

in  a  form  (such  as  a  binary  deck)  which  is  inconvenient  for  transfer 
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between  different  Fortran  and  machine  systems,  or  have  a  combination  of 
the  above  problems.   It  seemed  worthwhile,  therefore,  to  write  a  program 
usable  for  both  large  and  small  angular  momentum  values  in  Fortran  so 
that  it  may  be  used  on  essentially  any  reasonably  large  modern  computer 
by  at  most  a  few  changes  in  input/output  statements,  dimension  statements, 
and  function  definitions. 

2.   Description  of  the  Mathematical  Problem 

The  methods  for  calculating  the  3j  and  6j  coefficients  are  similar. 
The  9j  coefficient  is  calculated  as  a  sum  of  products  of  6j  coefficients. 
We  begin  with  a  discussion  of  the  calculation  of  the  6j  coefficients. 

A  general  formula  for  the  Wigner  6j  coefficient,  expressing  it  as  a 
single  sum,  was  first  given  by  Racah  [l]  and  is  given  as  Eq.  (6.3.7), 
page  99  in  Edmonds  [2]: 

x  v  r (-i)Z  (*+i>! 

ZiL(a-Ji-Ja-js)!    (z-Jx-Xs-4)!    (z-J&i-  js-4  )  • 
z 

x  1 

(z-ix-ig-js)!      (jl+j2+Xi+4|-z)!      (J3+J3+^S+X3-Z)I 

x  — — - — \ — —  I  (1) 

(a+b-c)!(a-b+c)!(-a+b+c)!  i 


where  A(abc)   = 


■J 


(a+b+c+l)i 

The  sum  is  over  all  non-negative  integer  values  of  z  such  that  no  facto- 
rial has  a  negative  argument.   The  function  A(abc)  is  defined  to  be  zero 
for  any  values  of  its  arguments  for  which  any  of  the  factorials  appearing 
therein  would  have  a  negative  argument.   Physically,  a  non-zero  value  of 

A  corresponds  to  two  angular  momentum  vectors  being  of  such  length  that 
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it  is  possible  for  them  to  couple  to  (or  form  a  triangle  with)  the  third. 
Any  combination  of  integer  and  ha  If- integer  angular  momentum  vectors 
which  would  lead  to  factorials  of  half-integers  in  the  formulae  is  also 
ruled  out.   Since  the  expression  on  the  right  hand  side  of  Eq.  (1)  is 
manifestly  symmetric  with  respect  to  any  permutation  of  the  indices  1, 
2  and  3,  there  is  no  necessity  or  gain  in  rearranging  the  coefficients 
using  its  symmetry  properties  (Ref.  2,  pp.  92-95,  and,  in  particular,  Eqs . 
(6.2.4)  and  (6.2.5)  on  pp.  94,95)  before  calculating  the  coefficient. 
Asymmetric  formulae  due  to  Racah  [l]  were  considered  for  the  calculation 
but  were  found  to  involve  series  with  the  same  number  of  terms  as  in  Eq. 
(1)  but  with  a  smaller  variation  in  the  size  of  the  terms,  which  makes  the 
errors  due  to  cancellations  between  successive  terms  in  the  alternating 
series  larger  (discussed  below).   An  additional  disadvantage  of  the  asym- 
metric formulae  is  that,  in  order  to  minimize  the  number  of  terms  in  the 
series,  the  rows  and  columns  of  the  symbol  must  be  rearranged  before  cal- 
culation, which  increases  the  computer  time  required. 

The  most  convenient  way  of  handling  the  large  numbers  represented  by 
the  factorials  is  to  take  the  logarithms  of  the  factorials.   Since  their 
calculation  is  relatively  time  consuming,  they  are  calculated  and  stored 
only  once  for  each  factorial,  the  first  time  called.   The  successive 
terms  in  the  series  in  Eq.  (1)  alternate  in  sign  and  vary  rather  slowly 
near  the  maximum  value,  leading  to  the  possibility  of  cancellation  errors 
due  to  addition  and  subtraction  of  nearly  equal  numbers.   Typical  behavior 
for  the  absolute  magnitudes  of  successive  terms  is  shown  in  Fig.  1.   Since 
the  range  of  the  absolute  values  of  the  terms  in  the  series  may  be  of  the 
order  of  1080  (e.g.,  for  j  and  I   of  the  order  of  60)  and  the  terms  near 


the  maximum  value  of  z,  z    ,  are  much  larger  for  large  angular  momenta 

max 

than  those  near  the  minimum  value  of  z,  z  .  ,  we  sum  the  series  downward 

nun 

starting  from  z   ,  and  cut  off  the  series  when  we  approach  terms  near 
b  max  r 

z  J      which  are  negligible  (see  Fig.  1).   The  first  term  in  the  series 
min  o  o  o 

(for  which  z  =  z    )  is  normalized  to  10~   .   This  allows  an  increase  by 

max  ■* 

CO 

a  factor  of  more  than  10   in  the  size  of  individual  terms  without  over- 

—  OR 

flow  and  permits  the  series  to  be  cut  off  at  10    without  problems  of 
underflow  in  double  precision.   The  normalized  series  is  then  summed 
directly  in  double  precision  arithmetic  and  multiplied  by  the  normaliz- 
ing coefficient,  which  is  the  product  of  the  value  of  the  term  for  which 

z  =  z    and  the  A's,  both  of  which  are  calculated  (in  version  I  in  single 

max  ° 

precision,  in  version  II  in  double  precision)  using  the  logarithms  of  the 

factorials  stored  in  the  subroutine  S6J: 

,Ji    Ja    js-i 

Ui    4    XsJ   =   ^hUU)   MJxJtaia)    A(4ij84)    Adiiaja) 

(-l)Zraax  (z   +1)' 
max 


(z   -ai)J  (z   -as)I  (z   -a3)I  (z   -a4)! 
max        max        max        max 


(bi-z   )!  (ba-z   )I  (b3-z    Ni 
max        max        max)T 


z 
max 

z    .    +1 
min 

X 

k=l 

\ 

(1') 


As  used  here,  single  precision  implies  8  significant  decimal  digits 
accuracy,  double  precision  means  16  significant  digits. 


where  ax   =   ji+js+j3 

*2     =    Jl+^2+4 
*3     =    ^l+J2+-«3 

a*   =  I-l+Zs+U 

bs     =    J2+J3+4+^3 

ba    =   ja+j  i+^3+^i 

z  =  rain    (bi  ,  bp  ,  b«  f 

max  c   A      s '   a  J 

Zmin  =  maX    fal  »a2 '^  ,a*^ 


ux   =    10 


,-15 


'(z        +l-a!-k)(z        +l-ag-k)(z        +l-a3-k)(z        +l-a4-k) 
max  max  max  max 


(z        +2-k)(b1-z        +k)(bo-z        +k)(ba-z        +k) 
L     max  max  max  max 


To  calculate  9j  coefficients  we  use  a  sum  of  products  of  6j  co- 
efficients, Eq.  (6.4.3)  on  p.  101  of  Ref.  1,  written  here  with  a  slight 
change  in  notation: 

■  -  2  (-l)aK  (2K+1) 

K 


jl  J2 

Jl2 

J3  J4 

J34 

J13 J24 

J 

,-Jl        J3     J1S'|      rJ2        J4     J24i      rjl2     J34     J      i 
tj24     J        *       J     tjs        H       J34J     U  ji        j2J 


(2) 


The  "coupling"  or  "triangle"  conditions  require  that  the  j's  in 
each  row  and  that  each  column  in  the  9j  symbol  form  a  triangle.   The  sum 
over  H  extends  from  the  largest  of  |ji-j|,  |j3-J24|»  an<^  |J2~J34|  to  tne 
smallest  of  ji+j,  J3+J24,  and  J2+J34.   If  the  "triangle"  conditions  are 
satisfied  for  the  9j  coefficients,  the  "triangle"  conditions  are  automa- 
tically satisfied  for  the  6j  coefficients  in  the  series  in  Eq.  (2). 


If  we  take  the  9j  coefficient  to  be  defined  as  the  sum  of  products  of 
three  6j  coefficients  with  the  parameters  distributed  as  in  (2),  then  the 
number  of  terms  in  the  sum  will  clearly  be  given  by 

N  =  H      -  H  .    -1 

max    min 
where  k^^  =  least  of  |ja+J34,  ji+j»  J3+js4| 

and  Hmin  =  greatest  of  jjjs- J34  | ,  |  ji- j  ! ,  |j3"ja4!} 

However,  using  the  symmetry  properties  (Ref.  4)  a  given  9j  symbol  may  be 
transformed  into  one  which  differs  in  value  by  at  most  a  sign,  but  in 
which  the  same  elements  are  differently  arranged,  and  having  a  different 
set  of  6j  coefficients  in  the  definition  in  (2),  and  hence  in  general  a 
different  number  of  terms,  N.   Of  the  72  permutations  given  by  the  symme- 
try properties  only  six  give  distinct  sets  of  three  6j  coefficients  in  (2); 
for  each  of  the  others  the  6j  coefficients  are  found  to  be  identical  to 
those  of  one  of  these  six,  utilizing  the  symmetry  relations  of  the  6j 
symbol  (invariance  with  respect  to  permutations  of  the  columns  and  in- 
variance  also  with  respect  to  simultaneous  inversion  of  any  two  columns). 
Using  the  notation  of  Ref.  4  we  may  choose  for  the  six  coefficients, 

1  2    3-n     r5    2    8^     rl   7   4- 
456  ><  417  W  396 
789'     U    3   9  '     U    8   5 

and 

2  1  3  -n  /-3  2  K  ^13  2 
546  i  654W465 
8   7   9  J  U   8   1J     W   9   8- 

of  which  the  first  three  are  identically  equal,  and  the  next  three  are 


R  9 

equal  to  the  first  multiplied  by  (-1)  ,  where  R  =  £  j  .  .   Again  in  the 

notation  of  Ref .  4,  the  first  is  in  the  class  Ci ,  to  which  there  corres- 
ponds the  unit  operator,  the  second  and  third  are  in  the  class  Qg  and 
correspond  to  the  permutations  (15)(38)(67)  and  (27X34M59)  respectively. 
The  last  three  are  in  the  class  C7   and  correspond  to  the  permutations 
(12M45M78),  (13)(46)(79),  and  (23)  (56M89) ,  respectively. 

For  speed  of  computation  we  want  the  series  for  the  9j  coefficient 
to  have  as  few  terms  as  possible  and  thus  we  compute  the  number  of  terms 
in  the  series  for  each  of  the  six  9j  coefficients  given  above  and  perform 
the  actual  computation  with  that  series  for  which  N  is  smallest.   It 
should  be  noted  that  this  procedure  may  result  in  a  much  shorter  series 
than  that  obtained  following  Ref.  2  (p.  106)  where  it  is  stated  that  the 
symbol  should  be  arranged  so  that  the  smallest  argument,  j  .  ,  does  not 
fall  in  the  positions  3,  5  or  7,  in  which  case  the  sum  over  K  has  at  most 

2i  ,   +1  terms.   Thus,  for  example,  the  coefficient 
Jmin 


17  11  12 
50  40  10 
65  50   15 


is  already  in  the  form 


such  that  the  smallest  argument,  j   =  10  is  not  in  the  position  3,  5  or 

6 

7,  and  in  this  form  the  series  (2)  has  20  terms  (k  .   =  2,  H    =  21). 
'  min       max 

However,  the  permutation  (15) (38) (67)  does  not  change  the  value  of  the 
9j  symbol  but  takes  it  into  the  form 


40  11  50 
50  17  65 
10   12   15 


,  for  which  the  series  (2)  has 


only  2  terms  (x    =  54,  K    =  55),  and  in  this  symbol  the  smallest 
J  min        max 


'  11    17 

12" 

40  50 

10 

.50  65 

15. 

argument  is_  in  the  position  7.  We  may  also  obtain  a  permutation  in  which 
the  smallest  argument  is  not  in  the  position  3,  5  or  7,  and  for  which  the 
series  also  has  only  2  terms:  viz.  (12) (45) (78) ,  which  takes  the  original 
9j  into  the  form 


without  change 


of  sign  since  R  =  i  j   is  even.   Thus  we  see  1)  that  the  procedure  used 

i=l 
in  the  current  program  may  reduce  the  number  of  terms  in  the  series  con- 
siderably and  2)  the  position  of  the  smallest  argument  in  the  symbol  is 
not  of  importance  in  determining  the  number  of  terms  in  the  series. 

Orthonormality  relations  for  testing  the  6j  coefficients  were  ob- 
tained from  Eq.  (6.2.9)  on  page  96  of  Ref.  2  (written  in  a  form  slightly 
more  convenient  for  our  purposes): 

hhhj  Ihhk  J    J3'J3 

and  for  the  9j  coefficients  from  Eq.  (6.4.6)  on  page  103  of  the  same 

reference : 

2      (2j12+l)(2j34+l)(2j13+l)(2j34+l) 

Jl2  J34 


J13  j  13  J24  J  24 


In  the  summation,  j12  runs  over  all  integral  (or  half- integral)  values 
to  which  ji  and  js  can  couple;  j34  runs  over  all  values  to  which  j3  and 
j4  can  couple. 
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If  the  Racah  W-coef f icient  is  desired,  it  may  be  obtained  from  the 
relation  (Eq.  (6.2.13),  p.  97  of  Ref.  2): 

w(jij2/s^;j3-e3)  =  (-D        -jjj  g  gl       (5) 

For  the  Wigner  3j  coefficient,  a  relation  combining  Eqs.  (3.6.11) 
and  (3.7.3)  on  pages  45  and  46  of  reference  2  is  used: 


Ji  Js  J3\  Jl-J3_m3 

jni  mg  1B3  /    m^  +mg  +m3  , 0 


((ji+js-ja)!  (j-i--ja+ja)i  (-Ji+js+Ja)J   I 

X    L  (J1+J3+J3+DJ  J 

X  I  (jj+mx  )•  C  j i - mx  ) !  ( J3+m8  )J  ( jg-rag  ) !( J3+m3  )J  (j3-m3  )  J 

x  Z  (-l)z  ^ 

z  zI    (  Jl+J3-J3-z)Hj1-ml-z)!(  j3+ms-z)I(J3-j2+m1+z)! 

x  71 r TTT  <6) 

The  orthonormality  relation  for  the  3j  symbols  used  for  testing  the 
accuracy  of  the  programs  is  Eq.  (3.7.8)  p.  47  of  Ref.  2: 

/Jl     J3     -J3\-/Jl     J3     J3     \ 

(2J3+D   2  )  (  ,)  =  6.    .    ,  6  ,   6(J!JaJ3)        (7) 

mlm3      ^ml     Hfe     ^3  '      NIDi     mg     m3     /  J3  J3  11131113 


The  mathematical  problems  involved  in  calculating  3j  coefficients 
for  large  values  of  the  arguments  are  similar  to  those  discussed  above  in 
connection  with  the  6j  coefficient.   We  therefore  do  not  repeat  the  main 
discussion  here. 


3.   Description  of  the  Programs 

Summary.   The  programs  were  originally  written  in  Fortran  II  for 
use  under  Fortran  Monitor  on  an  IBM  7094  computer.   Fortran  IV  programs 
have  been  added  since  this  language  is  being  adopted  increasingly  due  to 
its  greater  flexibility  and  speed.   Fortran  II  programs  are  listed  in 
Appendix  I,  and  Fortran  IV  in  Appendix  II. 

The  complete  package  consists  of  four  Fortran  functions,  F3J,  F6J, 
F9J  and  S6J.  F3J  (for  the  3j  coefficients)  is  independent  of  the  other 
three  except  in  the  Fortran  IV  version,  where  it  shares  a  common  memory 
block  with  S6J  in  which  the  logarithms  of  the  large  numbers  are  stored. 

F3J  performs  the  angular  momentum  coupling  tests,  tests  to  see  if 
mi+mg+ma  =  0,  and  tests  whether  j   and  m  are  either  both  integral  or 
both  half-integral  for  i  =  1,  2,  3.   If  any  of  these  tests  is  failed, 
the  subprogram  returns  a  zero  to  the  calling  program.   If  not,  the  sub- 
program proceeds  with  the  calculation  of  the  coefficient. 

F6J  performs  the  "triangle"  tests  to  see  if  the  Wigner  6j  coefficient 
is  zero  due  to  angular  momentum  coupling.   If  the  coefficient  is  zero, 
the  subprogram  returns  a  zero  to  the  calling  program.   If  not,  F6J  calls 
subprogram  S6J  which  calculates  the  coefficient. 

F9J  performs  the  triangle  tests  for  the  Wigner  9j  coefficient,  re- 
turns a  zero  if  the  arguments  fail  the  triangle  tests,  or  calculates  the 
coefficient  by  using  the  sum  of  a  series  of  products  of  6j  coefficients, 
Eq.  (2),  which  are  calculated  by  subprogram  S6J.   The  internal  sub- 
program S6J,  which  in  general  should  not  be  used  independently  of  F6J 
or  F9J  (unless  the  calling  program  eliminates  all  coefficients  which 
are  zero  due  to  angular  momentum  coupling),  calculates  and  stores  the 
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logarithms  of  the  first  322  factorials  the  first  time  called  and  then 
uses  these  to  calculate  the  6j  coefficient  using  Eq.  (l'). 

Comments  on  the  two  versions  of  the  program.   Version  I  is  in  single 
precision  except  for  the  summation  of  the  alternating  series  in  Eq.  (l') 
in  S6J  and  Eq.  (6)  in  F3J,  for  which  double-precision  is  used  because  of 
strong  cancellations  in  this  series,  as  is  discussed  below.   Version  II 

is  internally  completely  in  double  precision  (S6J  returns  a  double  pre- 
cision number  to  either  F6J  or  F9J)  but  both  F6J  and  F9J  return  a  single 
precision  number  to  the  calling  program.   The  errors  involved  in  each  of 
these  versions  are  given  in  Table  2. 

Function  S6J.   Since  in  general  a  computer  program  will  require  both 
6j  and  9j  coefficients,  we  use,  for  economy  of  space  in  the  computer  mem- 
ory, a  single  program,  S6J,  to  calculate  the  6j  coefficients  and  to  store 
the  logarithms  of  the  factorials  required  for  calculating  both  the  6j  and 
9j  coefficients,  rather  than  carry  out  these  operations  in  two  places  - 
once  when  computing  the  6j  coefficients  and  again  when  computing  9j  co- 
efficients.  The  subprogram  F6J,  which  provides  the  "triangle"  tests  and 
then  calls  S6J,  may  be  omitted  if  the  calling  program  provides  these  tests, 
but  should  be  used  otherwise.   F9J  also  performs  the  triangle  tests,  thus 
calling  only  6j  coefficients  which  are  not  zero  due  to  angular  momentum 
coupling.   F9J  therefore  calls  S6J  directly,  and  does  not  call  F6J. 

S6J  calculates  and  stores  the  logarithms  of  the  factorials,  FL(N)  = 
log  (N-l)I,  the  first  time  called.   Since  S6J  may  be  called  by  F6J  or  F9J 
or  by  appropriate  calculation  programs,  it  is  necessary  to  determine 
within  the  subprogram  S6J  itself  whether  or  not  this  subprogram  has  been 
called  previously.   This  is  carried  out  by  testing  whether  or  not  the 
ad  hoc  variable  NCALL  is  equal  to  -1867,  a  negative  prime  number.   Since 
it  is  extremely  unlikely  that  the  value  in  the  memory  cell  corresponding 

to  NCALL  is  initially  equal  to  this  particular  value,  we  have,  to  begin 
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with,  NCALL  ^   -1867,  in  which  case  the  program  then  sets  NCALL  =  -1867 

and  proceeds  to  the  calculation  and  storage  of  log  (N-D!  for  N=l  to 

322.   On  any  subsequent  calling  of  the  subprogram  S6J  the  variable  NCALL 

is  therefore  found  to  be  equal  to  -1867,  in  which  case  the  calculation 

of  the  quantities  log  (N-D!  (computed  the  first  time  S6J  was  called)  is 

bypassed. 

The  D0   loop  ending   in  statement   30  defines   z  (the   smallest  value 

of   z  in  Eq.    (1))    to  be   equal   to   the    largest  ojf  the   four  quantities 

(Ji+J2+J3)>    (ji+^+^s),    (^i+J2+^3)   and   Ux+Xa+ja).      The  D0    loop  ending 

in  statement  51  defines  z    (the  largest  value  of  z  in  Eq.  (1))  to  be 

max  °  n 

equal   to  the   smallest  of   the   three   quantities    ( J1+J2 +^2+^3  ) »    ( J2+J3+^s+^3  ) 
and   ( js  +  j  1  +^3  +h  )  • 


aWe  need,  in  the  calculation  of  the  3 j ,  6j  and  9j  coefficients,  log  (nl) 

for  n  =  0,1,  ...  N,  a  total  of  N+l  logarithms.   For  the  3j  coefficient 

\sax   ms    nfe/,   N  =    1+Ji+J2+J3«      Thus   for  all  j      =  80,   we  need   242    loga-< 

rithms   in  F3J.      For  the   6j   coefficient  j  j1    •£    S f  the    largest  N  arises 

in  the   factor   (z        +1)1    in  Eq.    (l').      Thus   in  S6J,N  =    1+z  where  z 

max  max        max 

=  smallest  of  the  three  numbers  ji+J2+^i+^  »  J2  +  J3  +4s  +-^3  ,  J3+ji+^3+^i,  so 
that  for  j  =  80  we  need  322  logarithms.   In  F9J  the  extreme  9j  coefficient, 
in  which  all  angular  momenta  take  on  the  value  80,  requires  for  its  com- 
putation 6j  coefficients  of  the  form  j    «n  80 1  w^ere  °  =  K  =  1^°*   How~ 

ever,  for  all  these  6j  coefficients  z    =  320,  independent  of  k,  so 

J  max 

that  for  the  calculation  of  the  9j  coefficients  we  also  require  only  322 

logarithms. 

In  Fortran  IV,  these  logarithms  are  stored  in  common,  so  that  we 

calculate  322  logarithms  in  either  F3J  or  S6J,  whichever  is  called  first. 
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The  summation  of  the  series  in  (l')  is  cut  off  when  the  absolute 
value  of  a  term  is  less  than  10~35(10~iO  times  the  first  term).  This  is 
particularly  useful  for  the  calculation  of  coefficients  involving  large 
angular  momenta  and  having  many  terms  in  the  series  (l')»  in  which  case 
there  are  many  terms  of  absolute  value  less  than  10  ,  the  calculation 
of  which  would  only  be  a  waste  of  time  since  they  are  completely  negli- 
gible compared  to  the  sum  of  the  series  (see  Table  1). 

Further,  within  the  subprogram  S6J  there  appears  the  factor  P,  de- 
fined two  lines  before  statement  72  (PUDG  =  XnP),  which  is  the  factor 
appearing  in  front  of  the  sum  in  Bq.  (l')  (apart  from  the  term  10   ). 
For  sufficiently  large  angular  momenta,  P  may  be  very  small  (<  10   ) 
although  the  6j  coefficient,  S6J  =  P*S  (line  after  statement  75)  is  of 

16  x 

order  10"1  or  10"3 ,  since  S  (S  =  sum  of  the  series  X10    is  very  large 
(see  Table  1).   In  order  to  avoid  possible  underflow  (which  occurs  at 

-66  .8 

about  10~29  =  e      when  using  double  precision  arithmetic  and  at  about 
10"38  =  e"87*6  in  single  precision)  in  the  calculation  of  P  (P  =  EXP(PL0G) 
in  statement  75)  we  introduce  the  statement  before  line  72: 

IF  (PL0G  +  64)  72,  75,  75 
and  work  with  Q 

Q  =  PL0G  +64. 
Thus  we  avoid  underflow  in  P  provided 

inP  >  -  130.8  (P>  2  X  10"57). 
This  limit  appears  to  be  satisfactory;  the  smallest  value  of  P  that  we 
have  encountered  appears  in  the  calculation  of  the  6j  coefficient 

{so  m  to}'  £or  whlch  inP  =  -113-2- 
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From  the  statements  above,  Q  <  0,  and  from  the  statements  following  state- 
ment 72,  viz., 

Q  =  EXP(Q) 
S6J  =  Q*S, 
we  will  not  have  overflow  provided  S<10   .   This  condition  seems  to  be 
satisfied  if  the  angular  momenta  occurring  in  the  6j  or  9j  coefficients 
are  all  =  80.   The  largest  value  of  S  we  have  encountered  is  in  the  6j 


:ore 


coefficient   J 80  8  8  1  for  which  S  =  1.8X1034.   The  statement  bef< 
^80  80  80J 

line  73,  viz.,  IF(ABS(S6J)-1. )  73,  73,  74,  causes  the  actual  6j  coefficient 

—  ft  A 

(which  is  e   X  S6J  appearing  immediately  above)  to  be  set  equal  to  zero 
if  it  is  less  than  e~s4(  =  1.6X10""28 ) . 

As  is  shown  by  the  values  of  | sum/maximum  term|  in  Table  1,  there 
are  strong  cancellations  of  the  terms  in  the  alternating  series,  Eq.  (l')» 
for  the  6j  coefficients,  and  the  necessity  for  the  double-precision  arith- 
metic  with  accuracy  10    used  in  summing  the  series,  both  in  versions  I 
and  II,  is  clearly  indicated  for  large  values  of  the  angular  momenta. 
The  limitation  of  the  programs  to  angular  momenta  of  80  is  because  the 
I  sum/maximum  term |  ratio  is  approaching  the  uncertainty  of  the  double  pre- 
cision arithmetic  (see  Table  1).   To  calculate  correctly  coefficients  with 
individual  angular  momenta  greater  than  80,  triple  precision  arithmetic 
would  be  required.   It  should  be  noted  (see  Table  1),  in  this  connection, 
that  for  angular  momenta  ^  80,  the  6j  coefficients  are  in  general  of  the 
order  10"   or  larger  (except  for  accidental  zeros). 

The  use  of  dummy  variables  (e.g.,  JD1)  in  the  calling  sequence  and 

setting  of  the  true  variables  of  the  subprogram  equal  to  them  (e.g., 

Jl  =  JDl)  minimizes  the  space  occupied  by  the  subprogram  and  maximizes 
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the  speed  of  execution  by  avoiding  a  large  number  of  location  instruc- 
tions, one  for  each  time  the  variable  appears  in  the  subprogram. 

Function  F9J.   This  subprogram  performs  angular  momentum  coupling 
tests  on  the  rows  and  columns  of  the  9j  coefficient  (all  statements 
through  statement  30)  and  rearranges  the  arguments  of  the  9j  coefficient 
to  select  the  shortest  series  possible  for  Eq.  (2)  as  described  earlier. 
It  then  calculates  the  9j  coefficient  using  Eq.  (2)  with  the  6j  coefficients 
provided  by  S6J. 

Function  F3J.   This  subprogram  performs  angular  momentum  coupling 
tests,  then  determines  the  minimum  value  of  z  to  be  the  absolute  value 
of  the  most  negative  of  0,  js-jg+ini*  and  J3_ji-ma.   The  maximum  value  of 
z  is  the  smallest  of  Ji+J2_J3»  ji_mi>  &n<3  J2+nfe.   The  coefficient  is  then 
calculated  using  Eq.  (6). 

Calling  Instructions.   The  arguments  of  the  functions  are  integers 
equal  to  twice  the  value  of  the  angular  momentum  involved.   Examples: 

ANSWER  =  F6J(J1,J2,J3,L1,L2,L3) 

B  =  F9J(J1,J2,J12,J3,J4,J34,J13,J24,J). 
To  calculate  the  6j  coefficient 


{ 


L      1      1-} 

,   ?  l\   use  C0EFF  =  F6J(8, 7,5,6, 7,3) . 

->   2   2  J 


Loading.   To  calculate  3j  coefficients  only,  load  F3J. 

To  calculate  6j  coefficients  only,  load  F6J  and  S6J. 

To  calculate  9j  coefficients  only,  load  F9J  and  S6J. 

To  calculate  both  6j  and  9j  coefficients,  load  F6J, 
F9J  and  S6J. 

To  calculate  3 j ,  6j,  and  9j  coefficients  load  F3J,  F6J, 
F9J,  and  S6J. 

Notes.   Programs  to  test  the  subprograms  F9J,  F6J,  and  F3J,  either 

by  calculating  individual  coefficients  or  doing  an  orthonormality  check, 

or  both,  are  given  in  the  appendices.   These  programs  are  labelled  TST9J, 
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TST6J,  and  TST3J.   In  the  testing  programs,  the  input  and  output  values 
of  the  arguments  are  the  angular  momenta  in  floating  point  form  (e.g., 
•j  =  3.5).   The  first  data  card  for  these  programs  is  for  labeling  runs 
and  columns  1-72  may  be  used  for  this  purpose. 

Examples  of  how  to  call  the  subroutines  are  given  in  the  testing 
programs,  TST9J,  TST6J,  and  TST3J. 

Version  II  takes  about  50%  more  time  for  a  given  calculation  than 
version  I.   Times  are  of  the  order  of  0.1  to  1  sec  per  9j  and  0.001  to  0.01 
sec  per  6j  but  vary  considerably  depending  on  the  configuration  assumed 
and  values  of  the  angular  momenta. 

4.   Results  and  Conclusions 

The  memory  occupied  by  and  the  accuracy  of  the  subprograms  described 
here  are  summarized  in  Table  2.   The  errors  quoted  were  determined  by  com- 
parison with  exact  6j  coefficients  or  by  use  of  the  orthonormality  rela- 
tions, Eqs.  (3),  (4),  and  (7). 

It  is  concluded  that  the  programs  described  here  operate  satisfac- 
torily for  angular  momenta  of  80  or  less,  and  that  programs  not  using 
integer  or  double  precision  arithmetic  to  avoid  cancellation  errors  will 
probably  not  operate  satisfactorily  for  angular  momentum  values  above 
about  40. 
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Fig.  1.   Semi- logarithmic  plot  of  the  absolute  values  of  successive  terms 
in  the  series  of  Eq.  (l')  for  j  =  20  for  all  arguments.   Note 
that  the  z  =  60  and  z  =  61  terms  are  cut  off  because  they  are 


smaller  than  10 


-26 


ia 


Table  1 

Properties  of  the  series  in  Eq.  (l')i 

(1)  Magnitude  of  the  term  having  the  largest  absolute  value. 

(2)  The  absolute  value  of  the  ratio  of  the  series  sum  to  the 
term  in  (1)  above. 


(3)  The  value  of  the  6j  coefficient  for  which  all  six  arguments 
are  equal: 

J  J  Jl 

•J  J  y' 


{: 


The  equal  argument  case  has  the  largest  maximum  terms  and  the 
strongest  cancellations  in  the  series  (smallest  value  of  | sum/ 
maximum  term | ) . 

(1)  (2) 

vsum 


J 

1 
10 
20 
40 
60 
80 


maximum  term 


10 


-16 


2  X  10 


-11 


1  X  10 

6  X  10€ 

5  X  10 

5  X  10 


.-6 
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,30 


maximum  term 


8  X  10 


,-i 


3  X  10 


,-2 


4  x  10 


,-3 


4  x  io' 


-8 


4  x  10 


,-8 


4  X  10 


,-j-a 


(3) 
6j 


.1666667  x  KT 


.2919187  x  10 


-a 


.5029406  X  10 


,-2 


.1828307  X  10 
.1006635  X  10 
.6568373  X  10' 


,-2 


-2 


-3 
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Table  2 

Comparison  between  Version  I  and  Version  II  (Fortran  II) 
memory  requirements  and  errors  in  Fortran  II  programs 


Program 


F9J 

F6J 

S6J 

F3J 

Average  j 

~  8 

(oi 

•thonormality) 

~  20 

ii 

~  40 

ii 

~  80 

ti 

Average 

j 

~  8 

(coefficient) 

~  20 

(or 

■thonormality) 

~  40 

ii 

~  80 

n 

Memory  Requirements  (words) 


Version  I 

Version  II 

962 

1066 

332 

339 

1031 

1554 

993 

1462 

Error  in  91 

Coefficient 

Version  I 

Version  II 

1.5  X  10"B 

3  X  10"7 

8.5  X  10"e 

1  X  10"6 

2  X  10' 

3  X  10' 


,-* 


1  X  10" 


,-• * 


1  X  10 


,-s 


Error  in  6j  Coefficient 
Version  I       Version  II 


1  x  10 


,-6 


1  X  10 


i-6 


1  X  10 


,-5 


8  X  10 


,-6 


1  x  10 

1  X  10 

1  x  10 

4  x  10' 


,-8 


,-V 


,-7 


,-7 


Average  j 


8  (coefficient) 

20  (orthonormality) 
40        " 
80        " 


Error  in  3J  Coefficient 
Version  I       Version  II 


3  X  10 


,-7 


1  X  10" 


4  X  10 


.-4 


8  X  10 


r4 


2  X  10 
1  X  10" 

3  X  10' 
5  X  10" 


,-8 


r7 


,-7 
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Appendix  I.   Listing  of  Fortran  II  Programs 
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Listing  of  Version  I  Fortran  II  Programs 
LABEL 


CF9J 


FUNCTIUM  F9J(JD1,JD2,JD3,JD4,JD5,JD6,JD7,JD8,JD9) 
;      F9J  VERSION  I   CALLS  S6J   FORTRAN  II 

DIMENSION  MTRIAI 18 ) ,KN( 6) ,KX(6) , NN( 6) 

J1=J01 

J2=JD2 

J3=J03 

J4=J04 

J5=JD5 

J6=JD6 

J7=JD7 

J8=JD8 

J9=JD9 
:      ANGULAR  MOMENTUM  COUPLING  TESTS  FOR  9J  COEFFICIENT 

I=-J1+J2+J3 

11=1/2 

IF(I-2*I1)  1000,1010,1000 
101C  MTRIA( 1  )  =  I  1 

I=+J1-J2+J3 

11=1/2 

IF  I  1-2*11)  1000,1020,1000 
1020  MTRIA(2)=I1 

I=+Ji+J2-J3 

11=1/2 

IF (  1-2*1 1)  1000,  1030, 1000 
103C  MTRIA(3)=I1 

I=-Jt+J5+J6 

11=1/2 

IF (1-2*11)  1000,1040,1000 
1040  MTRIA(4)=I1 

I=+J4-JS+J6 

11=1/2 

IF(  1-2*11)  1000,1050,1000 
1050  MTRIA{5)=I1 

I=+J4+J5-J6 

11=1/2 

IF(  1-2*11)  1000,1360,1000. 
1060  MTRIA(6)=I1 

l=-J7+J8+J9 

11=1/2 

IF(  1-2*11)  1000,1070,1000 
1070  MTRIAI 7)=I 1 

I=+J7-J3+J9 

11=1/2 

IF(  1-2*11)  1000,1080,1000 
1080  MTRIA(6)=I1 

I=  +  J  7+J8-J9 

11=1/2 

IF  (1-2*11)  1000,1090,1000 
1090  MTRIA{9)=I1 

I=-J1+J4+J7 

I  1=1/2 

IF(  1-2*11)  1000,1100,1000 
1100  MTRIA(10)=I1 

I=+Jl-J4+J7 

11=1/2 

IF (  1-2*1 1 )  1000,1110,1000 
1110  MTRIA( 11 )=I1 

I=+J1+J4-J7 
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1120 


1130 


1140 


1150 


1160 


1170 


1180 


30 


35 


40 

50 


52 


11  =  1 
IF{  I 

MTRI 
I=-J 
11  =  1 
IF(  I 
MTRI 

I  =  +  J 

II  =  1 

I  F  (  I 
MTRI 
I=+J 

II  =  1 
IF(  I 
MTRI 
I=-J 
11=1 
IF(  I 
MTRI 
I=  +  J 
11  =  1 
IF(  I 
MTRI 
I=  +  J 
11  =  1 
IF(  I 
MTRI 
DO  3 
IF(M 
CONT 
KNU 
KN12 
KM  (3 
KN(4 
KN(5 
K  N  (  6 
KX(  1 
KX12 
KX(3 
KX(4 
KX(5 
KX(6 
DO  3 
NNIK 
KSIG 
I=XM 
DO  4 
IFII 
COMT 
KMIN 
KMAX 
GO  T 
J  =  J1 
J1  =  J 
J5  =  J 
J=J3 
J3=J 
J8  =  J 
J=J6 
J6  =  J 


8 


*I1)  1000,1120,1000 

12)=I1 

J5  +  J8 

*I1)  1000,1130,1000 

13)=I1 

J5  +  J8 

Ml)  1000,1140,1000 

14)=I1 

J5-J8 

*I1)  1000,1150,1000 
15  )  =  I  1 

J6+J9 

*I1)  1000,1160,1000 
16 )  =  1 1 
J6  +  J9 

*I1)  1000,1170,1000 

17  )  =  1 1 

J6-J9 


*I1)  10 
1 8  )  =  1 1 
N=lf18 

IA(N) ) 

U£ 

XMAX0F( 

XMAXOF ( 

XMAXOF { 

XMAXOF ( 

XMAXOF ( 

XMAXOF ( 

XMINOFl 

XMINOFl 

XMINOFl 

XMINOFl 

XMINOF ( 

XMINOFl 

K  =  l,6 

KX(K)-K 

1 

OFINNI 1 

K  =  l,6 

<1  (  K  )  )  40 

UE 

,N  I  K  )  + 1 

X(K)+1 

130,52, 


00, 1180,1000 


1000,30,30 


XA3SFU2 
XABSFl J2 
XABSFl J6 
XARSF ( Jl 
XABSFl J2 
XA3SFU3 
J2+J6, Jl 
J2+J7, J5 
J6+J7, J5 
J1+J6.J2 
J2+J4, J3 
J3+J5, Jl 

NIK) 

),NN(2), 

,50,40 


-J  6), XABSFl J  1- J 9), XABSFl J4-J8)  ) 
-J  7), XABSFl J  5- J 9), XABSFl J4-J3) ) 
-J7), XABSFl J5-J1  ), XABSFl J8-J3) ) 
-J  6)  , XABSFl J2-J9), XABSFl J  5- J  7)  ) 
-J  4)  , XABSFl  J 3-J7)  , XABSFl J6-J8)  ) 
-J  5), XABSFl  JI-J8), XABSFl J4-J9)  ) 
+J9, J4+J8) 
+J9, J4+J3) 
+J1, J8+J3) 
+J9, J5+J7) 
+J7, J6+J8) 
+J8, J4+J9) 


NN  (  3  )  t  NN  (  4  )  ,  AIN  (  5  )  f  NN  (  6  )  ) 


53,54,55,56)  ,K 
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53 


54 


55 


56 


120 


130 


200 


1000 
2  COO 

C 
C 
* 
CF6J 


J7  =  J 

GO  TO 
J  =  J2 
J2=J7 
J7=J 
J=J3 
J3=J4 
J4  =  J 
J=J5 
J5  =  J9 
J9  =  J 
GC  TU 
J=J1 
J1=J2 
J2=J 
J=J4 
J4  =  J5 
J5=J 
J  =  J7 
J7=J8 
J8=J 
GO  TO 
J  =  J1 
J1  =  J3 
J3  =  J 
J=J4 
J4  =  J6 
J6  =  J 
J  =  J7 
J7  =  J9 
J9  =  J 
GO  TO 
J  =  J2 
J2  =  J3 
J3=J 
J=J5 
J5  =  J6 
J6  =  J 
J  =  J8 
J8  =  J9 
J9  =  J 
KSIGN 
SUMMA 
SUM  =  0 
SIG=( 
FLK=K 
DO  20 
TERM  = 
1*S6J( 
FLK=F 
SUM=S 
F9J  =  S 
GO  TO 
F9J  =  0 
RETUR 
END 


13C 


13£ 


120 


120 


(1-X 
TION 
.0 

-1)** 
M 1 N 
0  K=K 
FLK*S 
J3»  J6 
LK+2. 
UM+TE 
UM*SI 

2  000 
»0 
M 


MODFI Jl+J2+J3+J4+J5+J6+J7+J8+J9,4) ) 
OF  SERIES  OF  EQUATION  (2) 

(KMIN-1 )*KSIGN 

M IN, KM AX, 2 

&J(Jl,J4,J7rJ8,J9,K-l)*S6J(J2,J5,J8,J4,K-l,J6) 
,J9,K-1, Jl, J2> 

0 

RM 

G 


LABEL 
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FUNCTION  F6J(JD1,J02,JD3,LD1,LD2,LD3) 
;      F6J  VERSION  I  CALLS  S6J  VERSION  I   FORTRAN  II 

DIMENSION  MED(1?) 

Jl=JOl 

J2=JD2 

J3=JD3 

L1=LD1 

L2=LD2 

L3=LH3 

ANGULAR  MOMENTUM  COUPLING  TESTS  FOR  6J  COEFFICIENT 

I=-J1+J2+J3 

11=1/2 

IF  (1-2*11)  1000,1010,1000 
1000  F6J=0. 

GO  TO  100 
1010  MED(1)=I1 

I=+J1-J2+J3 

11=1/2 

IF  (1-2*11)  1000,1020,1000 
1020  MED(2)=I1 

I=+J1+J2-J3 

11=1/2 

IF  (1-2*11)  1000,1030,1000 
1030  ME0(3)=I1 

I=-J1+L2+L3 

11=1/2 

IF  (1-2*11)  1000,1040,1000 
1040  MED(4)=I1 

I=+J1-L2+L3 

11=1/2 

IF  (1-2*11)  1000,1050,1000 
1050  MED(5)=I1 

I=+J1+L2-L3 

11=1/2 

IF  (1-2*11)  1000,1060,1000 
1060  ME0(6)  =  U 

I=-L1+J2+L3 

11=1/2 

IF  (1-2*11)  1000,1070,1000 
1070  MED(7)=I1 

I=+L1-J2+L3 

11=1/2 

IF  (1-2*11)  1000,1080,1000 
1080  MED(8)=I1 

I=+L1+J2-L3 

11=1/2 

IF  (1-2*11)  1000,1090,1000 
1090  MED(9)=I1 

I=-L1+L2+J3 

11=1/2 

IF    (1-2*11)     1000,1100,1000 
HOC    MED(10)  =  I1 

I=+L1-L2+J3 

11=1/2 

IF  (1-2*11)  1000,1110,1000 
1110  MED(11)=I1 

I=+L1+L2-J3 

11=1/2 

IF  (1-2*11)  1000,1120,1000 
1120  MED(12)=I1 
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10 
100 

c 
c 

• 
CS6J 

c 


50 
15 


20 

30 


40 
51 


DO  10  N=l,12 

IF  (MED(N))  1000,10,10 

CONTINUE 

F6J  =  S6JU1,J2,J3,L1,L2,L3) 

RETURN 
END 


LABEL 


WHETHER  TO  CALCULATE 
5,15,5 


FL(N)    S 


FUNCTION    S6J( JD1, JD2, JD3,LD1, LD2,LD3) 

VERSION    I       FORTRAN     II 

DIMENSION    MA(4),MB(3),MED(12),FL<322) 

J1=JD1 

J2=JD2 

J3=JD3 

L1=LD1 

L2=LD2 

L3=LD3 

DETERMINE 

IFINCALL+1867) 

NCALL=-1867 

CALCULATE    FL(N)    S 

FL( 1)=0.0 

FL(2)=0.0 

DO    50    N=    3,322 

FN=N-1 

FL(N)=FL(N-1)+L0GF(FN) 

MED(l)=(-Jl+J2+J3>/2 

MED(2)=(+Jl-j2+J3)/2 

ME0(3)=(+Jl+J2-J3)/2 

*IEDU)  =  (-Jl  +  L2+L3)/2 

MED(5)=<+Jl-L2+L3)/2 

MED(6)  =  UJl  +  L2-L3)/2 

MED{  H  =  (-Ll+J2+L3)/2 

MED(8)=(+Ll-J2+L3)/2 

MEQ(9)=(+Ll+J2-L3)/2 

ME0(10)=(-Ll+L2+J3)/2 

MED( 11  )=(+Ll-L2+J3)/2 

MED( 12  >=(+Ll+L2-J3)/2 

MA(1)=MED(1 )+MED(2)+MED<3) 

MA ( 2 ) =MED ( 4 ) +MED ( 5 )  +  MED  (  6  ) 

MA  J  3  >=MED(7)+MED(8 )+MED(9) 

MA(4)=MED( 1C )+MED( 11)+ME0( 12) 

MB(1  )=MA( 1)+MED(  12) 

MB(2)=MA( 1)+MED(4) 

MB(3)=MA(  1)+MED(  8  ) 

DETERMINE  MAXIMUM  OF  (J1+J2+J3) 

MAX  =  i*A(  1  ) 

DO  30  IM=2,4 

IF  (MAX-MA(N))  20,30,30 

MAX=MA(N) 

CONTINUE 

DETERMINE  MINIMUM  OF  ( J  1  + J2+L  1+L2  )  , 

MIN  =  MB(  1  ) 

DO  51  N=2,3 

IF  (MIN-MB(N))  51,51,40 

MIN=MB(N) 

CONTINUE 

KMAX=MIN-MAX 
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(Jl+L2+L3),(L1+J2+L3),(L1+L2+J3) 


( J2+J3+L2+L3) , ( J3+J1+L3+L1) 


55 


HI 
MI 
MI 
MI 
MI 
MI 
MI 
MI 
MI 
SU 
UK 
S  = 
IF 
DO 
UK 


60 
65 


70 


1 


72 


73 

74 
75 

78 

80 

90 


CU 
IF 
S  = 
S  = 

:a 

DE 
DO 
NU 
DE 
NU 
NU 
NU 
NU 
DE 
DE 
UL 

PL 
IF 
Q= 
Q= 
S6 
IF 
S6 
30 
S6 
GO 
P= 
S6 
MI 
IF 
S6 
CO 
RE 
EN 


NP1=MIN+1 

N1=MINP1-MA( 1 ) 

N2=MINP1-MA< 2) 

N3=MINP1-MA( 3) 

N4=MINP1-MA( 4) 

N5=MINP1+1 

N6  =  Ml>(1  )-MIN 

N7=MB(2)-MIN 

N8=MB(3)-MIN 

M  SERIES  IN  DOUBLE  PRECISION 

=1„E-15 
l.Lt-15 

(KMAX)  65,65,55 

60  K=1,KMAX 
=-UK*       FLOATF(MIN1-K)*FLOATF(MIN2-K)*FLOATF(MIN3-K)* 

FL0ATF(MIN4-K)/(FL0ATF(MIN5-K)*FL0ATF(MIN6+K)*FL0ATF(MIN7+K)* 

FL0aTF(MIN8+K ) ) 
T  OFF  SERIES  AT  l.OE-25 
(  A3SF(UK)-l.E-25)  65,60,60 
S+UK 

S*1.0E  +  1.5 

LCULATE  DELTA  FUNCTIONS 
LOG=0.0 

70  N=l,12 
M=MED(N) 

L0G=DEL0G+FL(NUM+1) 
M1=MA( l)+2 
M2=MA(2)+2 
M3=MA(3)+2 
M4=MA(4)+2 

L0G=D£L0G-FL(NUM1  )-FL  (NUM2  )-FL  (  NU/^3  )-FL  (  NUM4  ) 
L0G=0„5*0EL0G 
0G=FL(MIN5 )-FL(MINl)-FL(MIN2)-FL(MIN3)-FL(MIN4J-FL(MIN6+l)- 

FL(MIN7+1)-FL<MIN8+1 ) 
OG=0€LDG+ULOG 
(PLOG+64. )     72,75,75 
PLiJG  +  64, 
EXPF(  j) 
J  =  Q*S 
(Ai3SF(S6J)-l.)  73,73,74 

J  =  0  >  a 

TO  90 
J=S6J*EXPF(-64C ) 

TU  78 
EXPF(PLOG) 
J  =P*S 
N2=MIN/2 

(MIN-2*MIN2)  80,90,80 
J=-S6J 
NTINUE 
TURN 
D 


CF3J 


LABEL 

FUNCTION  F3J( JD1, JD2, JD3,MD1,MD2,MD3) 
F3J  VERSION  I  FORTRAN  II 
DIMENSION  MTRI  (9 ) ,FL( 242) 
J1=JD1 
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50 
15 


J2=JI>2 

J3  =  j  :>3 

M1=MD1 

M2=M02 

M3=MD3 

IF(NCALL+1867)5, 15, 5 

NCALL=-1867 

FL(1)=0. 

FL(2)=G. 

CO    5  0    N=3,242 

FN=N-1 

FL (N)=FLIN-1 )+LOGF(FN) 

I=J1+J2-J3 


11=1/2 

IFII-2*U)  1000,  101J,10C0 
1010  M,TRI(1)  =  I1 

I=J1-J2+J3 

11=1/2 

IF  <  1-2*  1 1 J  1000,1020,1000 
1020  MTRI(2)=11 

I=-JH-J?  +  J3 

11=1/2 

I  F  {  1-2*1 1)  1000,1030,1000 
1030  MTRI(3)=I1 

IF(M1+M2+M3)  100C,1040tlGOO 
1040  I=J1+M1 

11=1/2 

IF  (1-2*11)  1000,1050,1000 

MTRI (4  )  =  I1 

MTRI (5 ) =( Jl-Ml )/2 

I  =  J2  +  |V2 

11=1/2 

IF (1-2*11)  1000,1060,1000 

MTRI  (6)=I  I 

MTRI ( 7)=( J2-M2)/2 

I=J3+M3 

11=1/2 

I F ( 1-2* I  1  )  1000,1070,1000 

MTRI (8)=Il 

MTRI (9)=( J3-M3J/2 

DO  30  M=l,9 

IF(MTRKM))  1000,30,30 

CONTINUE 

IF  (J3-J2+M1)  41,45,45 

KMIN=-J3+J2-M1 

GO  TG  60 

KMIN=0 

IF(-J3+J1+M2  -KMIN)  80,80,70 

KMIN=-J3+J1+M2 

KMIN=  KMIN/2 

IFU2-J3+M1)     90,100,100 

KMAX=J1+J2-J3 

GO    TO    11C 

KMAX=J1-M1 

IF  (  J2  +  M2-KMAX)     120,130,130 

KMAX=J2+M2 

KMAX=KMAX/2 

M  IN1  =  MTRI I  1 )-KMIN+l 

MIN2=MTRI (5)-KMIN+l 

MIN3=MTRI (6J-KMIN+1 
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1050 
1060 

1070 

30 

40 

45 
60 
70 

80 

90 

100 
110 
120 
130 


MIN4=(  J3-J2  +  MD/2+KMIN 
MIM5=( J3-Jl-M2)/2+KMIN 
D  UK    =1.E-10 

D  S=1.E-10 

NCUT=0 

KMAX=KMAX-KMIN 
IFJKMAX)  165,  165,155 
155  00  160  K=1,KMAX 
D      UK=-UK*FL0ATF(MIN1-K)*FL0ATF(MIN2-K)»FL0ATF(MIN3-K)/ 
1  (FLOATF  (KMIN+K  )*FL0ATF(  MIN4+K  )  *F|_0ATF  (  M  IN5  +  K  )  ) 
IFUbSF<UK)-l.E30)  158,157,157 
D  157  UK=1.E-10*UK 
D      S=1.E-10*S 
NCUT=NCUT+1 
158  IF(ABSF(UK)-l.E-20)  165,160,160 
D  160  S=S+UK 
165  DEL0G=0. 

DO  170  N=l,9 
NUM=MTRI (N) 

170  DEL0G=DEL0G+FL(NUM+1) 
NUM=(Jl+J2+J3)/2+2 
DEL0G=G.5*(DEL0G-FL(NUM) ) 

UL0G=-FL(KMIN+1)-FL(MIN1)-FL<MIN2)-FL(MIN3)-FL(MIN4+1)-FL(MIN5+1) 
PLOG=DELOG+ULOG 

IF(PLOG+80.)  172,171,171 

171  IF(NCUT)  175,175,172 

172  SIGN=SIGNF(1.,S) 
S=ABSF(S) 

SL0G=L0GF(S)+FL0ATF(NCUT+1)*L0GF( l.F+10) 
F3J=SIGN*EXPF(SL0G+PL0G) 

30  T0178 
175  S=S*1.E+10 

P=EXPF(PL0G) 

F3J=P*S 
178    NUM=KMIN+( Jl-J2-M3)/2 

IF(XMODF(NUM  ,2))  180,190,180 
180  F3J=-F3J 
190  CONTINUE 

GO  TO  2  000 
1000  F3J=08 
2000  RETURN 

END 
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Listing  of  Version  II  Fortran  II  Programs 
LABEL 


CF9J 


FUNCTION  F9J( JD1,JD2,JQ3, JD4,JD5,JD6, JD7,JD8,JD9) 
I      F9J  VERSION  II  CALLS  S6J   FORTRAN  II 

DIMENSION  MTRIA( 18) ,KNl 6) ,KX(6) ,NN{ 6) 

J1=JD1 

J2=J02 

J3=JD3 

J4=JD4 

J5=JD5 

J6=JD6 

J7=JD7 

J8=JD8 

J9=J09 
;      ANGULAR  MOMENTUM  COUPLING  TESTS  FOR  9J  COEFFICIENT 

I=-J1+J2+J3 

11=1/2 

IF (  1-2*1 1 )  1000, IOIOtIOOO 
1010  MTRIA(1)=I1 

I=+J1-J2+J3 

11=1/2 

I F (  1-2* 1 1 J  1000,1020,1000 
1020  MTRIA(2)=I1 

I=+J1+J2-J3 

11=1/2 

IF(  1-2*1 1 )  1000,1030,1000 
1030  MTRIA(3)=I1 

I=-J4+J5+J6 

11=1/2 

IF (  1-2*1 1)  1000,1040,1000 
1040  MTRIA(4)=I1 

I=+J4-J5+J6 

11=1/2 

I F {  1-2* I  1 )  1000,1050,1000 
1050  MTRIA(5)=I1 

I=+J4+J5-J6 

11=1/2 

IF (  1-2*1 1)  1000,1060,1000 
1060  MTRIA(6)=I1 

I=-J7+J8+J9 

11=1/2 

IF (  1-2*1 1)  1000,1070,1000 
1070  MTRIA(7)=I1 

1=+J7-J8+J9 

11=1/2 

IF (  1-2*1 1 )  1000,1080,1000 
1080  MTRIA(8)=I1 

I=+J7+J8-J9 

11=1/2 

IF (  1-2*11 )  1C00, 1090,1000 
1090  MTRIA(9)=I1 

I=-J1+J4+J7 

11=1/2 

IF (  1-2* I  1  )  1000,1100,1000 
1100  MTRIA(10)=II 

I=+J1-J4+J7 

11=1/2 

IF (  1-2*1 1  )  1000,1110,1000 
1110  MTRIA(11)=I1 

I=+J1+J4-J7 

30 


1120 


1130 


1140 


1150 


1160 


1170 


118C 


30 


35 


40 
50 


52 


II 
IF 
MT 
1  = 
II 
IF 
MT 
1  = 
II 
IF 
MT 
1  = 
II 
IF 
MT 
1  = 
II 
IF 
MT 
1  = 
II 
IF 
MT 
1  = 
II 
IF 
MT 
DO 
IF 
CO 
KM 
KN 
KN 
KN 
KN 
KN 
KX 
KX 
KX 
KX 
KX 
KX 
DO 
NN 
KS 
1  = 
DO 
IF 
CO 
KM 
KM 
SO 
J  = 
Jl 
J5 
J  = 
J3 
J8 
J= 
J6 


=  1/2 

(  1-2*11 

»IA( 12) 

-J2+J5+ 

=  1/2 

(  1-2*11 

R  I A  (  1 3  ) 

+J2-J5+ 

=  1/2 

(  1-2*11 

RIA( 14) 

+J2+J5- 

=  1/2 

(  1-2*11 

RIA{ 15) 

-J3+J6+ 

=  1/2 

(  1-2*11 

RIA(16) 

+J3-J6+ 

=  1/2 

(  1-2*11 

RIA(17) 

+J3+J6- 

=  1/2 

(1-2*11 

RIA( 18) 

30  N  =  l 
(MTRIA( 
NTINUE 
(1 )=XMA 
(2)=XMA 
(3)=XMA 
(4)=XMA 
(  5  )  =  X  A  A 
(6)=XMA 
<1)=XMI 
(2)=XMI 
(3)=XMI 
(4)=XMI 
(5)=XMI 
(6)=XMI 

35  K  =  l 
(K)=KX( 
IGN=1 
XMINGFl 

40  K=l 
(  I-NN(K 
NTINUE 
IN=KN(K 
AX=KX(K 

TO (130 
Jl 
=  J5 
=J 
J3 
=  J8 
=  J 
J6 
=  J7 


)  1000,1120,1000 

=  11 

J8 

)  1000,1130,1000 
=  11 

J8 

)  1000,1140,1000 

=  11 

J8 

)  1000,1150,1000 

=  11 

J9 

)  1000,1160,1000 

=  11 

J9 

)  1000,1170,1000 

=  11 

J9 

)  1000,1180,1000 

=  11 

,18 

N))  1010,30,30 


X0F(XARS 
XOFUABS 
XOFUABS 
XOFUABS 
X0F(XA8S 
XOFIXARS 
NOF ( J2+J 
NOF( J2+J 
NOF( J6+J 
NOF( Jl+J 
NCF( J2+J 
NO  F.I  J  3+ J 
,6 
K)-KN(K) 

NN( 1),NN 

,6 

)  )40,5% 


F( J2-J6) ,  XABSF { J1-J9),XABSF( J4-J8)  ) 
F( J2-J7),XA8SF( J  5- J9  ) , AABSF ( J4-J3)  ) 
F( J6-J7) , XA3SF (J5-J1  )  , XABSF ( J8-J3 )  ) 
F( Jl-J6)» XABSF (J2-J9)  , XABSF ( J 5- J 7  )  ) 
F( J2-J4 ),XABSF(J3-J7),XaBSF(J6-J8)) 
F( J3-J5) , XABSF ( J1-J8), XABSF ( J4-J9)  ) 
6, J1+J9, J4+J8) 
7, J5+J9, J4+J3) 
7, J5+J1, J8+J3) 
6, J2+J9, J5+J7) 
4, J3+J7, J6+J8) 
5, J1+J8,J4+J9) 


(  2  )  ,  NN  (  3  )  ,  NN  (  4  ) ,  NN  (  5  ) ,  :mN  {  6  )  ) 
40 


)  +  l 

)  +  l 
,52,53,54,55,56) ,K 
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13! 


120 


J7  =  J 

GO  TO  130 

53  J  =  J2 
J2=J7 
J7=J 
J  =  J3 
J3=J4 
J4  =  J 
J=J5 
J5  =  J9 
J9  =  J 
GO  TO 

54  J=J1 
J1=J2 
J2=J 
J=  J4 
J4=J3 
J5=J 
J=J7 
J7  =  Jd 
J8=J 
GO  TO 

55  J=J1 
J1=J3 
J3  =  J 
J=J4 
J4  =  Jo 
J6  =  J 
J=J7 
J7  =  J9 
J9  =  J 
GO  TO 

56  J=J2 
J2=J3 
J3  =  J 
J=  J5 
J5=J6 
J6  =  J 
J=J8 
J8  =  J9 
J9=J 

120  KSIGN=( l-XMODF( J 1+J 2+J3+J4+J5+J 6+J7  +  J 8+ J9,  4 )  ) 
C      SUMMATION  OF  SERIES  OF  EQUATION  (2) 
D  130  SUM=0.C 


l?; 


SIG=(-1)**(KMIN-1)*KSIGN 

Y  =  Go 

Y=KMIN 

PLK  =  Y 

DO    200    K  =  KMIN,MAX,2 

TERM=FLK*S6J(Jl,J4,J7„J8tJ9,K-l)*S6J<J2,J5,J8tJ4,K-l,J6) 
1»S6J( J3, J6, J9,K-1, Jl, J2 ) 

FLK=FLK+2.0 
200  SUM=SUM+TERM 

SUM=SUM*SIG 

F9J=SUM 

GO  TO  2000 
1000  F9J=^,0 
2000  RETURN 

END 
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c 

CF6J 

C 


1000 
1010 

1020 

1030 

1040 
1050 
1060 
1070 
1080 
1090 
1100 
1110 


LABEL 

FUNCTI 
F6J  V 
DIMEMS 
J1=JD1 
J2=JD2 
J3=JU3 
L1=L01 
L2=L02 
L3=LD3 
ANGULA 
I=-J1+ 
11=1/2 
IF  (  I- 
F6J=Lo 
GC  TO 
MED{1) 
I=+J1- 
11=1/2 
IF  (  I- 
MED(2) 
I=+J1+ 
11=1/2 
IF  (  I- 
MED(3) 
I =-J 1+ 
11=1/2 
IF  (  I- 
MED(4) 
I=+J1- 
11=1/2 
IF  (  I- 
ME0(5) 
I=+J1+ 
11=1/2 
IF  (  I- 
MEO(b) 
I=-L1+ 
11=1/2 
IF  (  I- 
MEO( 7) 
I=+L1- 
11=1/2 
IF  (I- 
MED(B) 
I=+L1+ 
11=1/2 
IF  (  I- 
MED19) 
I=-L1+ 
11=1/2 
IF  (I- 
MEO( 10 
I=+L1- 
11=1/2 
IF  (  I- 
MEDdl 
I=+L1+ 


ON  F6JUDI  ,  JD2,JD3,LD1,LD2,LD3) 

ERSION  II  FORTRAN  II  CALLS  S6J  VERSION  II 

ION  ME0(12) 


R  MOMENTUM  COUPLING  TESTS  FOR  6J  COEFFICIENT 
J2  +  J3 

2*11)  1000,1010,1000 

100 

=  11 

J2  +  J3 

2*11)  1000,1020,1000 

=  11 

J2-J3 

2*11)  1000,1030,1000 
=  11 
L2  +  L3 

2*11)  1000,1040,1000 

=  11 
L2+L3 

2*11)  1000, 1G50, 1000 

=  11 
L2-L3 

2*11)  ICO"',  1060, 1000 
=  11 
J2  +  L3 

2*11)  1000,1070,1000 

=  11 
J2  +  L3 

2*11)  1000,1080,1000 

=  11 

J2-L3 

2*11)  1000,1090,1000 

=  11 
L2+J3 

2*11)  1000,1100,1000 

)=I1 

L2  +  J3 

2*11)  1000,1110,1000 

)  =  I1 

L2-J3 
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1120 


1C 


100 

c 

c 

* 
CS6J 

c 

D 


50 
15 


WHETHER  TO  CALCULATE  FL(N)  S 
,15,5 


20 

30 


11=1/2 

IF  (I-2»I1)  1000,1120,1000 

MED(12)=I1 

DO  10  N=l,12 

IF  (MED(N))  1000,10,10 

CONTINUE 

A6J=S6J(J1,J2,J3,L1,L2,L3) 

F6J=A6J 

RETURN 

END 


LABEL 

FUNCTION  S6J( JD1 , JD2 , JD3 , LD1 , LD2, LD3 ) 

S6J  VERSION  II  FORTRAN  II 

DIMENSION    FL(322) 

DIMENSION  MA(4),MB(3),MED(  12) 

J1=J01 

J2=J02 

J3=JD3 

L1=LD1 

L2=L02 

L3=LD3 

DETERMINE 

IF(NCALL+1867) 

NCALL=-1867 

CALCULATE  FL(N)  S 

FN  =  0. 

FL(1)=0«0 

FL(2)=0oQ 

DO  50  N=  3,322 

A=N-1 

FN  =  A 

FL(N)=FL(N-1)+L0GF( FN) 

MED(l)=(-Jl+J2+J3)/2 

ME0(2)=(+Jl-j2+J3)/2 

MED(3)=(+Jl+J2-J3)/2 

MED(4)=(-Jl+L2+L3)/2 

ME0(5)=(+Jl-L2+L3)/2 

MED(6)=(+Jl+L2-L3)/2 

MED(7)=(-Ll+J2+L3)/2 

MED(8)=(+Ll-j2+L3)/2 

MED(9)=(+Ll+J2-L3)/2 

MED(10)=(-Ll+L2+J3)/2 

MED( 11 )={+Ll-L2+J3)/2 

MED(12)=(+Ll+L2-J3)/2 

MA( 1 )=MED( 1)+MED(2)+MED(3) 

MA(2)=MEDU)+ME0(5)+MED(6) 

MA(3)=MED(7)+ME0(8)+MED{9) 

MA(4)=MED( 10)+MED( 11 )+MED( 12) 

MB(1)=MA( 1 )+MED( 12) 

MB(2)=MA( 1 )+MED( 4) 

MB(3)=MA( 1 )+MED( 8) 

DETERMINE  MAXIMUM  OF  ( J  1  + J2  + J3)  ,  ( Jl +L2  +  L3 ) 

MAX  =  MA(  1) 

DO  30  N=2»4 

IF  (MAX-MA(N))  20,30,30 

MAX=MA( N) 

CONTINUE 

34 


(L1+J2+L3), (L1+L2+J3) 


C      DETERMINE  MINIMUM  OF  ( J 1+J2+L 1+L2) ,  ( J2+J3+L2+L3) , ( J3+J1+L3+L1 ) 

MIN=MB( 1 ) 

DO  51  N=2,3 

IF  (MIN-MBlN)J  5  1,51,40 
40  MIN  =  MB(ISI) 
51  CONTINUE 

KMAX=MIN-MAX 

MINP1=MIN+1 

MIN1=MINP1-MA{ 1) 

MIN2=MINP1-MA<2) 

MIN3=MINP1-MA{  3) 

MIN4=MINP1-MA<4) 

MIN5=MINP1+1 

MIN6  =  M6(  D-MIN 

MIN7=MB(2)-MIN 

MIN8=MB(3)-MIN 
C      SUM  SERIES  IN  DOUBLE  PRECISION 
D      UK=loE-15 
D      S=1.0E-15 

IF  (KMAX)  65,65,55 
55  DO  63  K=1,KMAX 
D      UK=-UK*       FL0ATF(MIN1-K)*FL0ATF(MIN2-K)«FL0ATF(MIN3-K)* 

1  FLCUTF(MIN4-K)/(FL0ATF(MIN5-K)*FL0ATF(MIN6+K)*FL0ATF(MIN7+K)* 

2  FL0ATF(MIN8+K  )  ) 

C      CUT  OFF  SERIES  AT  l.OE-25 

IFlABSF(UK)-l.E-25)  65,60,60 
D   60  S=S+UK 
D   65  S=S»1.0E+15 
C      CALCULATE  DELTAS 
D      DELOG=0.0 

DO  70  N=l,12 

NUM=MED(N) 
D   70  DEL0G=DEL0G+FL(NUM+1) 

NUM1=MA( l)+2 

NUM2=MA(2)+2 

NUM3=MA(3)+2 

NUM4=MA(4)+2 
D      DEL0G=DEL0G-FL(NUM1)-FL(NUM2)-FL(NUM3)-FL(NUM4) 
D      DELOG=0,5*DELOG 
D      UL0G=FL(MIN5)-FL(MIN1)-FL(MIN2)-FL(MIN3)-FL(MIN4)-FL(MIN6+1)- 

1    FL(MIN7+1)-FL(MIN8+1) 
D      PLOG=DELOG+ULOG 

IF(PL0G+64.)  72,75,75 


D 

72 

Q=PL0G  +  64. 

D 

Q=EXPF(0) 

D 

S6J=Q*S 

IF(A6SF(S6J)-1.)      73,73,74 

0 

73 

S6J=0. 
GO  TO  90 

D 

74 

S6J=S6J»EXPF(-64o ) 
GO  TO  78 

D 

75 

P=EXPF(PLOG) 

D 

S6J  =P*S 

78 

MIN2=MlM/2 

IF  (MIN-2*MIN2)  80,90,80 

D 

80 

S6J=-S6J 

90 

CONTINUE 

D 

RETURN 
END 
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c 

CF3J 

C 
D 


D 

0       50 
15 


1010 

1020 

1030 
1040 

1050 
1060 

1070 

30 

40 

45 
60 
70 
80 


LABEL 

FUNCTI 

F3J    VE 

DIMENS 

D I  MENS 

J1=JD1 

J2=JD2 

J3=J03 

M1=MD1 

M2=MD2 

M3=MD3 

DETERM 

IF(NCA 

NCALL= 

FL<1)= 

FL(2)= 

DO    50 

A=N-1 

FN=A 

FL(N)= 

I=J1+J 

11=1/2 

IF ( 1-2 

MTRK1 

I=J1-J 

11=1/2 

IF{ 1-2 

MTRK2 

I=-J1+ 

11=1/2 

I F  i 1-2 

MTRI  (3 

IFJM1+ 

I=Jl+M 

11=1/2 

1 F { 1-2 

MTRI (4 

MTRI (5 

I=J2+M 

11=1/2 

IF(  1-2 

MTRK6 

MTRI (7 

I=J3+M 

11=1/2 

IF( 1-2 

MTRI (8 

MTRI (9 

DO    30 

IF(MTR 

CONTIN 

IF     ( J3 

KMIN=- 

GO    TO 

KMIN=0 

IF(-J3 

KMIN=- 

KMIN=K 


ON    F3J( JD1 ,JD2, JD3,MD1,MD2,MD3) 
RSION     II     FORTRAN     II 
ION    FH242) 
ION    MTRK9) 


INE    WHETHER    TO    CALCULATE    FL ( N )     S 

LL+1867J5, 15,5 

-1867 

0. 

0. 

N=3,242 


FL(N-1)+L0GF(FN) 
2-J3 

♦  id   iooo,  mio,iooc 

)  =  I1 
2  +  J3 

*I1)  1000,1020,1000 

)=I1 

J2  +  J3 

*I1)  1000,1030,1000 

)  =  U 

M2+M3)  1000,1040,1000 

1 

♦III  1000,1050,1000 

)  =  I1 

)=(Jl-Ml)/2 

2 

*I1)  1000,1060,1000 

)  =  I1 

)=( J2-M2)/2 

3 

*I1)  1000,1070,1000 

)  =  I1 

)=( J3-M3J/2 

N=l,9 

I(N))  1000,30,30 

UE 

-J2+M1)  40,45,45 

J3+J2-M1 

60 

+J1+M2  -KMIN)  80,80,70 
J3+J1+M2 

MIN/2 
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IFU2-J3  +  M1)     90,100,100 
90    KMAX=J1+J2-J3 
30    TO     11C 
100    KMAX=J1-M1 

110  IF( J2+M2-KMAX)  120,130,130 
120  KMAX=J2+M2 
13C  KMAX=KMAX/2 

MIN1=MTRI ( 1 J-KMIN+1 
MIN2=MTRI ( 51-KMIN+l 
MIN3=MTRI (6)-KMIN+l 
MIN4=( J3-J2+Ml)/2+KMIN 
MIN5=lJ3-Jl-M2)/2+KMIN 
C      SUM  SERIES  IN  DOUBLE  PRECISION 
D      UK  =1.E-10 
D      S=1.E-10 
NCUT=0 

KMAX=KMAX-KMIN 
IF(KMAX)  165,  165,155 
155  DO  160  K=1,KMAX 
D      UK=-UK*FL0ATF(MIN1-K)*FL0ATF(MIN2-K)*FL0ATF(MIN3-K)/ 
1(FL0ATF(KMIN+K)*FL0ATF(MIN4-+K)*FL0ATF(MIN5+K)  ) 
IF(ABSF(UK)-1.E30)  158,157,157 
D  157  JK=1.E-10*UK 
D      S=1.E-1C*S 
NCUT=NCUT+1 
158  IF(ABSF(UK)-l.E-20)  165,160,160 
D  160  S=S+UK 

c    calculate  delta  functions 

D  165  DELOG=0, 

D      DO  170  N=l,9 

NUM=MTRI (N) 
D  170  DEL0G=DEL0G+FL(NUM+1) 

NUM=( Jl+J2+J3)/2+2 
D      DELOG=0.5*(DELOG-FL(NUM)  ) 

D      UL0G=-FL(KMIN+1 )-FL ( MINI )-FL( MIN2 )-FL (MIN3 )-FL ( MIN4  + 1 )-F L( MIN5+1 ) 
D      PLOG=DELOG+ULOG 

IF(PL0G+80.)  172,171,171 
171  IF(NCUT)  175,175,172 
D  172  SIGN=SIGNF(1.,S) 
D      S=ABSF(S) 

D      SL0G=LUGFIS)+FL0ATF(NCUT+1)*L0GF(1.E+10) 
D      F3J=SIGN*EXPF(SL0G+PL0G) 

GO  T0178 
D  175  S=S*1.E+10 
D      P=EXPF(PLOG) 
D      F3J=P«S 

178  NUM=KMIN+( J1-J2-K3J/2 

IF(XMODF(NUM  ,2))  180,190,180 
D  180  F3J=-F3J 
190  CONTINUE 

GO  TO  2000 
1000  F3J=0. 
2000  RETURN 

END 
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Listing   of    Fortran    II   Testing   Programs 
#  LABEL 

CTST9J 
C  TEST    OF    F9J    FUNCTION    SUBPROGRAM       FORTRAN    II 

DIMENSION  AJ19) 
C      RUN  LABELLING, USE  FIRST  DATA  CARD  COLUMNS  1-72 
READ  INPUT  TAPE  5,10 
10  F0RMATC72H 

1  ) 

WRITE  OUTPUT  TAPE  6,10 
C      READ  MODE*  =1  DO  COEFFICIENTS  ONLY,  =2  DO  BOTH  COEFFICIENTS  AND 
C      ORTHONORMALITY,  =3  DO  ORTHONORMALITY  ONLY 
READ  INPUT  TAPE  5, 20, MODE 
IF(M0DE-2)  28,28,52 
C      READ  NMAX,NUMBER  OF  INDIVIDUAL  COEFFICIENTS  TO  BE  CALCULATED 
28  READ  INPUT  TAPE  5,20,NMAX 

WRITE  OUTPUT  TAPE  6,30,NMAX 
20  F0RMATI7I10) 

30  F0RMATI24H  NUMBER  OF  CALCULAT IONS=, 13 ) 
32  WRITE  OUTPUT  TAPE  6,35 

35  F0RMAT(80H       Jl       J2       J3       J4       J5       J6       J7 
1       J8       J9       9J) 
DO  50  N=1,NMAX 
C      READ  TRUE  Jl , J2, J3, J4, J  5, J6, J7, J8, J9  FLOATING  POINT  E.G.  3/2  =  1,5 
READ  INPUT  TAPE  5  ,40 , < A J ( I ) , 1  =  1 , 9) 
40  F0RMAT(9F8.1) 

J1=2.»AJ(1)+.001 
J2=2.*AJ(2)+.001 
J3=2.*AJ(3)+.001 
J4=2.*AJ<4)+.001 
J5=2.*AJ(5)+.001 
J6=2.*AJ(6)+.001 
J7=2„»AJ<7)+.001 
J8=2„*AJ(8)+.001 
J9=2.*AJ(9)+.001 

ANSWER=F9JJJl,J2,J3,J4,J5,J6,J7,jS,J9) 
WRITE  OUTPUT  TAPE  6,45, ( A JII ) , I =1 , 9 ) , ANSWER 
45  F0RMAT(9F8.1,E15.7) 
50  CONTINUE 

52  IF(M0DE-2)  200,55,55 
C      9J  ORTHONORMALITY  CHECK 
C      READ  NMAX,  NUMBER  OF  NORMALIZATION  RUNS  TO  BE  MADE 

55  READ  INPUT  TAPE  5, 20, NMAX 
WRITE  OUTPUT  TAPE  6, 30, NMAX 
DO  100  N=1,NMAX 

WRITE  OUTPUT  TAPE  6,56 

56  F0RMAT(72H       Jl       J2       J3       J4      J13    J13PR      J24 
1    J23PR        J  ) 

C      READ  TRUE  Jl , J2, J3, J4, J  13, J13PR, J24, J24PR, J    FLOATING  POINT 

READ  INPUT  TAPE  5, 40 ,TJ 1, TJ2, TJ 3,TJ4, TJ 13, TJ13PR, TJ24, TJ24PR,TJ 

WRITEOUTPUTTAPE  6,40,TJ 1,TJ2, TJ3, TJ4, TJ13.TJ13PR, T J24,TJ24PR,TJ 

J12MIN=2.*ABSF(TJ1-TJ2)+.001 

J12MAX=2.»!TJ1+TJ2)+.001 

J34MIN=2.«ABSF(TJ3-TJ4)+,001 

J34MAX=2.»(TJ3+TJ4)+.001 

J1=2.*TJ1+.001 

J2=2.»TJ2+.001 

J3=2.*TJ3+.001 

J4=2.*TJ4+.001 

J13=2.«TJ13+.001 

J13PR=2.0*TJ13PR+.001 
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J24=2.*TJ24+.001 
J24PR=2.0«TJ24PR+.001 
J=2.»TJ+.001 
C      9J  ORTHONORMALITY  CHECK  BY  EQUATION  (4) 
SUM=0. 

DO    100    J12=J12MIN,J12MAX,2 
DO    100    J34=J34MIN,J34MAX,2 
IFU12  +  J34-J)    100,70,70 
70    IF<XABSF( J12-J34)-J)80,80,100 
80    TERM=FLOATF{ J12+1 ) *FLOATF< J34+1 )*FLOATF( J13+1 )*FLOATF( J24+1 ) 

1*F9J(J1,J2,J12,J3,J4,J34,J13,J24,J)*F9J(J1,J2,J12,J3,J4,J34,J13PR, 
2J24PR,J) 
SUM=SUM+TERM 
TJ12=0.5*FLOATF( J12) 
TJ34=0.5*FL0ATF< J34) 

WRITE  OUTPUT  TAPE  6, 60, TJ 12,TJ34, TERM, SUM 
60  F0RMAT(2F6.1,2E15.7) 
100  CONTINUE 
200  CALL  EXIT 
END 
C 
C 

*      LABEL 
CTST6J 
C      PROGRAM  FOR  TESTING  F6J  FUNCTION  SUBPROGRAM      13  JULY  1965 

DIMENSION  AJ<6) 
C      RUN  LABELLING, USE  FIRST  DATA  CARD  COLUMNS  1-72 
READ  INPUT  TAPE  5,10 
10  F0RMAT172H 

1  ) 

WRITE  OUTPUT  TAPE  6,10 
C      READ  MODE,  =1  DO  COEFFICIENTS  ONLY,  =2  DO  BOTH  COEFFICIENTS  AND 
C      ORTHONORMALITY,  =3  DO  ORTHONORMALITY  ONLY 
READ  INPUT  TAPE  5, 20, MODE 
IF1M0DE-2)  28,28,52 
C      READ  NMAX, NUMBER  OF  INDIVIDUAL  COEFFICIENTS  TO  BE  CALCULATED 
28  READ  INPUT  TAPE  5, 20, NMAX 

WRITE  OUTPUT  TAPE  6, 30, NMAX 
20  F0RMAT(7I10) 

30  F0RMAT(24H  NUMBER  OF  CALCULAT IONS=, 13 ) 
32  WRITE  OUTPUT  TAPE  6,35 

35  F0RMAT(55H       Jl       J2       J3       LI       L2       L3      6J> 
DO  50  N=1,NMAX 
C      READ  TRUE  Jl , J2, J3 , LI , L2, L3  IN  FLOATING  POINT  FORM  E.G.  3/2=1.5 
READ  INPUT  TAPE  5,40, ( AJ ( I ) , 1=1 ,6) 
40  F0RMAT(9F8.1)   " 
J1=2.*AJ(1)+.001 
J2=2.*AJ<2)+.001 
J3=2.*AJ(3)+.001 
J4=2.*AJ(4)+.001 
J5=2.»AJ(5)+.001 
J6=2.*AJ(6)+.001 
ANSWER=F6J{ Jl, J2,J3,J4, J5,J6) 

WRITE  OUTPUT  TAPE  6 ,45, ( AJ( I ) , I =1,6 ) , ANSWER 
45  F0RMAT(6F8.1,E15.8) 
50  CONTINUE 

52  IF(M0DE-2)  200,55,55 
C      ORTHONORMALITY  CHECK 

C      READ  IMAX,  NUMBER  OF  ORTHOGONALITY  RUNS  TO  BE  MADE 
55  READ  INPUT  TAPE  5, 20, IMAX 

39 


WRITE  OUTPUT  TAPE  6,30, IMAX 
DO  150  1=1, IMAX 
WRITE  OUTPUT  TAPE  6,56 
56  F0RMAT(48H       Jl       J2       J3       LI       L2     J3PR) 
C      READ  TRUE  Jl , J2, J3, LI ,L2, J3PRIME  FLOATING  POINT  AS  ABOVE 
READ  INPUT  TAPE  5,40,TJ 1.TJ2, TJ3, TL 1, TL2,TJ3PR 
WRITEOUTPUTTAPE  6, 40, TJ 1,TJ2, TJ3.TL 1,TL2,TJ3PR 
J1=2.*TJ1+.001 
J2=2.*TJ2+.001 
J3=2.*TJ3+.001 
L1=2.*TL1+.001 
L2  =2.*TL2  +.001 
J3PR=2.*TJ3PR+.001 

IF(XABSF(J2    -LI    )-XABSF(Jl    -L2    >)  80,80,90 
80  NMIN  =  XABSFU1    -L2    )+l 

GO  TO  100 
90  NMIN=XABSFU2    -LI    )+l 
100  IFJJ2    +L1    -Jl    -L2    )  110,110,120 
110  NMAX=J2    +L1    +1 

GO  TO  130 
120  NMAX=J1    +L2    +1 
C      6J  ORTHONORMALITY  BY  EQUATION  (3) 
130  SUM=0. 

DO  150  N=NMIN,NMAX,2 

TERM=FLOATF(N)*FLOATF{ J3PR  +1 )*F6J ( Jl, J2, J3,L1 ,L2,N-1 ) *F6J( Jl , J2 , 
1J3PR,L1,L2,N-1) 
SUM=SUM+TERM 
140  F0RMAT(3H  N=,I3,6H  TERM=, E16.8, 5H  SL'M  =  ,E16.8) 

WRITE  OUTPUT  TAPE  6, 140 ,N, TERM, SUM 
150  CONTINUE 
200  CALL  EXIT 
END 
C 
C 

*      LABEL 
CTST3J 

C      PROGRAM  FOR  TESTING  F3J  FUNCTION  SUBPROGRAM   FORTRAN  II 
DIMENSION  AJ(3),AM(3) 
READ  INPUT  TAPE  5,10 
10  FORMAT (72H 

1  ) 

WRITE  OUTPUT  TAPE  6,10 
READ  INPUT  TAPE  5, 20, MODE 
WRITE  OUTPUT  TAPE  6, 25, MODE 
25  F0RMAT(6H  M0DE=,I1) 
I,F(M0DE-2)  28,28,52 
28  READ  INPUT  TAPE  5,20,NMAX 

WRITE  OUTPUT  TAPE  6,30,NMAX 
20  FORMAT(TIIO) 

30  F0RMAT(24H  NUMBER  OF  CALCULAT IONS=,  13 ) 
32  WRITE  OUTPUT  TAPE  6,35 

35  F0RMAT(55H       Jl       J2       J3       Ml       M2       M3      3J) 
DO  50  N=1,NMAX 

READ  INPUT  TAPE  5,40, { AJ( I ) , 1=1 , 3) , < AM( I ) , 1=1, 3) 
40  F0RMAT(9F8.1) 

J1=2.*AJ<1)+.001 
J2=2.*AJ(2)+,001 
J3=2.*AJ<3)+.001 
M1=2.001*AM{ 1) 
M2=2.001*AM(2) 
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M3=2.001*AM(3) 

ANSWER=F3Jt J1,J2,J3,M1,M2,M3) 

WRITE  OUTPUT  TAPE  6 ,45, ( A J ( I ) , 1=1, 3> , ( AM( I ) , 1=1 ,3 ) , ANSWER 
45  FORMAT(6F8.1,E15.8) 
50  CONTINUE 
52  IFIM0DE-2)  200,55,55 

ORTHONORMALITY  CHECK 

55  READ  INPUT  TAPE  5,20, IMAX 
WRITE  OUTPUT  TAPE  6,30,  IMAX 
DO  150  1=1, IMAX 

WRITE  OUTPUT  TAPE  6,56 

56  F0RMAT<48H       Jl       J2        J3    J3PR       M3     M3PR) 
READ  INPUT  TAPE  5,40 ,TJ 1 ,TJ2, TJ3,TJ3PR, TM3, TM3PR 
WRITEOUTPUTTAPE  6,40  ,TJ1, TJ2, TJ3,TJ3PR, TM3, TM3PR 
J1=2.0001*TJ1 

J2=2.0001*TJ2 

J3=2.0001*TJ3 

J3PR=2.0001*TJ3PR 

M3=2.0001»TM3 

M3PR=2.0001*TM3PR 

SUM=0. 

M1ST0P=2*J1+1 

M2STOP=2*J2+l 

DO  150  M=1,M1ST0P,2 

DO  150  N=1,M2ST0P,2 

M1=M-J1-1 

M2=N-J2-1 

TERM=FL0ATF<J3+1)*F3J<J1,J2,J3,M1,M2,M3)»F3J(J1,J2,J3PR,M1,M2,M3PR 
1) 

SUM=SUM+TERM 

IF(TERM-l.E-lO)  150,150,145 
140  F0RMAT(4H  M1=I3,4H  M2=,I3,6H  TERM=, E16. 8, 5H  SUM=,E16.8) 
145  WRITE  OUTPUT  TAPE  6, 140, Ml, M2,  TERM, SUM 
150  CONTINUE 
200  CALL  EXIT 

END 
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Appendix  II.   Listing  of  Fortran  IV  Programs 
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Listing   of   Version   I   Fortran   IV  Programs 


$IBFTC    F9J  LIST 

FUNCTION    F9J(JD1,JD2,JD3,JD4,JD5,JD6,JD7,J0  8,JD9) 
C  F9J    VERSION     I       CALLS    S6J       FORTRAN     IV 

DIMENSION    MTRIA( 18) ,KN( 6)  ,KX(6)  ,NN( 6) 

J1=JD1 

J2=J02 

J3=JU3 

J4=JD4 

J5=JD5 

J6=JD6 

J7=J07 

J8=JD8 

J9=JD9 
C  ANGULAR    MOMENTUM    COUPLING    TESTS    FOR    9J    COEFFICIENT 

I=-J1+J2+J3 

11=1/2 

IF  ( 1-2*11 )  1000,1010,1000 
1010  MTRIA(1)=I1 

I=+J1-J2+J3 

11=1/2 

IF  (  1-2*1 1 )  1000,1020,1000 
1020  MTRIA(2)=I1 

I=+J1+J2-J3 

11=1/2 

IF(I-2*I1)  1000,1030,1000 
1030  MTRIA(3)=I1 

I=-J4+J5+J6 

11=1/2 

I  F{  1-2* 1 1 )  1000,1040,1000 
1040  MTRIA(4)  =  U 

I=+J4-J5+J6 

11=1/2 

IF(I-2*I1)  1000,1050,1000 
1050  MTRIA(5)=I1 

I=+J4+J5-J6 

11=1/2 

IF(I-2*I1)  1000,1060,1000 
1060  MTRIA(6)  =  U 

I=-J7+J8+J9 

11=1/2 

IF(I-2*I1)  1000,1070,1000 
1070  MTRIA(7)=I1 

I=+J7-J«+J9 

11=1/2 

IF(I-2*I1)  1000,1080,1003 
1C80  MTRIA(8)=I1 

I=+J7+J8-J9 

11=1/2 

IFU-2*I1)  1000,1090,1000 
1090  MTRIA(9)=I1 

I=-J1+J4+J7 

11=1/2 

IF ( 1-2*11 )  1000,1100,1000 
1100  MTRIA(10)=I1 

I=+J1-J4+J7 

11=1/2 

IF(I-2*U)  1000,1110,1000 
1110  MTRIA(11)=I1 

I=+J1+J4-J7 

11=1/2 
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1120 


1130 


1140 


1150 


1160 


1170 


30 


I F (  1-2*  11} 

MTRIA(12)= 
I=-J2+J5+J 
11=1/2 
IF(  1-2*11) 
MTRIA( 13)= 
I=+J2-J5+J 
11=1/2 
IF(I-2*I1) 
MTRIA{ 14)= 
I=+J2+J5-J 
11=1/2 
IF(  1-2*11) 
MTRIA( 15)= 
I=-J3+J6+J 
11=1/2 
IF<  1-2*11) 
MTRIA( 16)= 
I=+J3-J6+J 
11=1/2 
IF (  1-2*1 1) 
MTRIA(17)= 
I=+J3+J6-J 
11=1/2 
IF(  1-2*11) 
1180  MTRIA(16)= 
DO  30  N=l, 
IF<MTRIA(N 
CONTINUE 
KN( 1)=MAX0 
KN(2)=MAX0 
KN(3)=MAX0 
KN<4)=MAX0 
KN(5)=MAX0 
KN(6)=MAX0 
KX( 1)=MIN0 
KX(2)  =  MI,M0 
KX(3)=MIN0 
KX(4)=MIN0 
KX<5)=MIN0 
KX(6)=MIN0 
DO  3  5  K=l 
NN(K)=KX(K 
KSIGN=1 
I=MIN0(NN{ 
DO  40  K=lf 
IF( I-NN(K) 
CONTINUE 
KMIN=KN(K) 
KMAX=KX (K) 
GO  T0(13C, 
J=J1 
J1=J5 
J5  =  J 
J=J3 
J3  =  J8 
J8  =  J 
J  =  J6 
J6  =  J7 
J7  =  J 


lOCOt 1120,1000 
II 
8 


1000, 
II 


1000, 

II 


10C0, 

II 

9 
1000, 

II 

9 
1000, 

II 

9 


1130,  100C 


1140,1000 


1150,1000 


1160, 1000 


1170,1000 


IABS( J6-J8) ) 
IABS( J4-J9) ) 


35 


40 
50 


52 


1000, 1180,1000 
II 
18 
))  1000,30,30 

( IABS( J2-J6),  IABS( J1-J3),  IABS( J4-J8) ) 

( IABSl J2-J7), IABS( J5-J9), IABS( J4-J3) ) 

( IABS( J6-J7),  IABS( J5-J1), IABS( J8-J3) ) 

( IABS< J6-J1), IABS( J2-J9),  IABS( J5-J7) ) 

{ IABSl J 2- J4),  IA8S( J3-J7) 

( lABSf J3-J5),  IABS( J1-J8) 

( J2+J6.J1+J9, JA+J8) 

( J2+J7, J5+J9, J4+J3) 

[ J6+J7, J5+J1, J8+J3) 

( J1+J6,J2+J9,J5+J7) 

( J2+J4,J3+J7, J6+J8) 

( J3+J5, J1+J8, J4+J9) 

6 

)-KNlK ) 

1) ,NN(2) ,NN(3),NN(4),MN(5),NN(6) ) 

6 

)40, 50,40 

+  1 
+  1 
52, 53, 54, 55, 56), K 
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GO  TO  13u 

53  J=J2 
J2=J7 
J7  =  J 
J=J3 
J3  =  J4 
J4  =  J 
J=J5 
J5=J9 
J9  =  J 

GO  TO  130 

54  J=J1 
Jl=J2 
J2  =  J 
J=J4 
J4  =  J5 
J5=J 
J=J7 
J7=J8 
J8=J 

GO  TO  120 

55  J=J1 
J1=J3 
J3  =  J 
J=J4 
J4  =  J6 
J6  =  J 

J  =  J7 

J7  =  J9 

J9  =  J 

GO  TO  120 

56  J=J2 
J2  =  J3 
J3  =  J 
J=J5 
J5  =  J6 
J6=J 
J=J8 
J8  =  J9 
J9=J 

120  KS IGN= ( l-MOD( J  1  + J 2  + J  3+ J 4+ J 5  + J 6+ J 7+ J 8  + J9, 4) ) 
C      SUMMATION  OF  SERIES  OF  EQUATION  (2) 
130  SUM=OoO 

SIG=(-1)**(KMIN-1)*KSIGN 
FLK=KMIN 

00  200  K=KMIN,KMAX,2 

TERM=FLK*S6J (J1,J4,J7TJ8,J9,K-1)*S6J(J2,J5,J8,J4,K-1,J6) 
l*S6J(J3,J6»J9,K-ltJl,J2) 
FLK=FLK+2.0 
200  SUM=SUM+TERM 
F9J=SUM«SIG 
GO  TO  2000 
1000  F9J=0.0 
2000  RETURN 
END 
C 
C 
SIBFTC  F6J      LIST 

FUNCTION  F6JIJD1 ,JD2,JD3,LD1,LD2,LD3) 
C      VERSION  I  F6J  FUNCTION  CALLS  S6J   FORTRAN  IV 
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1000 
1010 

1020 

1030 

1C40 

1050 

1060 

1070 

1080 

1090 

1100 

1110 

1120 


DIMENS 

J1=J01 

J2=JD2 

J3=JD3 

L1=LD1 

L2=LD2 

L3=L03 

ANGULA 

I=-J1+ 

11=1/2 

IF  (I- 

F6J=0. 

GO  TU 

MED(l) 

I=+J1- 

11=1/2 

IF  (  I- 

MED(2) 

I=+J1+ 

11=1/2 

IF  (  I- 

MED(3) 

I=-J1+ 

11=1/2 

IF  (I- 

MED(4) 

I=+J1- 

11=1/2 

IF  (  I- 

MED(5) 

I=+J1+ 

11=1/2 

IF  (  I- 

MED(6) 

I=-L1+ 

11=1/2 

IF  (  I- 

MED(7) 

I=+L1- 

11=1/2 

IF  (  I- 

med(«) 

I=+L1+ 
11=1/2 
IF  (  I- 
MED(9) 
I=-L1+ 
11=1/2 
IF  (  I- 
MED( 10 
I=+L1- 
11=1/2 
IF  (  I- 
VIEDI  11 
I=+L1+ 
11=1/2 
IF  (  I- 
MEO( 12 
DO  10 
IF  (ME 


ION  MED(12) 


R  MOMENTUM  COUPLING  TESTS  FOR  6J  COEFFICIEMT 
J2  +  J3 

2*1 1 )  1000,1010,1000 

0 

100 

=  11 

J2  +  J3 

2»I1)  1000,1020,1000 

=  11  • 

J2-J3 

2»I1)  1000,1030,1000 
=  11 
L2  +  L3 

2*11)  1000,1040,1000 
=  11 
L2  +  L3 

2*11)  1000,1050,1000 

=  11 

L2-L3 

2*11)  1000,1060,1000 
=  11 
J2  +  L3 

2*11)  1000,1070,1000 
=  11 
J2  +  L3 

2*11)  1000,1080,1000 

=  11 

J2-L3 

2*11)  1000,1090,1000 
=  11 
L2  +  J3 

2*11)  1000,1100,1000 
)  =  I1 
L2  +  J3 

2*11)  1000,1110,1000 

)  =  I1 

L2-J3 

2*11)  1000,1120,1000 

)  =  I1 

N=l,12 

0(N))     1000,10,10 
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10  CONTINU 

F6J=S6J 
100  RETURN 

END 
C 
C 
SIBFTC  S6J 

FUNCTIO 
C      VERSION 

DIMENSI 

COMMON/ 

DOUBLE 

DFLOAT( 

Jl=JDl 

J2  =  Jl;2 

J3=JD3 

L1=L01 

L2=LD2 

L3=LD3 
C      DETERMI 

IF(NCAL 

5  NCALL=- 

C      CALCULA 

FL(1)=U 

FL(2)=0 

DO  50  N 

FN=N-1 
50  FL(N)=F 
15  HED(1)= 

MED(2)= 

MED(3)= 

MED(4)= 

MED(5)= 

MED(6)= 

MED(7) = 

MED(8)= 

MED(9)= 

MED (10) 

MED(ll) 

MEDU2) 

MA(1)=M 

MA(2)=M 

MA(3)=M 

MA(4)=M 

MB(  1)=M 

MB(2)=M 

MB(3)=M 
C      DETERMI 

MAX=MA( 

DO  30  N 

IF  (MAX 

MAX=MA( 

CONTINU 

DE TERM IN 

MIN=M3( 1 

DO  51  N= 

IF  (MIN- 

MIN=M8(N 

CONTINUE 

KMAX=MIN-MAX 


( Jl , J2rJ3,Ll,L2,L3) 


LIST 
N  S6J( JD1 , JD2, JD3,LD1,LD2,LD3) 

I   FORTRAN  IV 
ON  MA(4) ,M3(3) ,MED( 12) 
FACT/FLI322) ,NCALL 
PRECISION   UK,S,OFLOAT 
I)  =  I 


20 
30 


40 
51 


NE  WHETHER  TO  CALCULATE  FL(N)  S 

L+1867)  5,15,5 

1867 

TE  FL(N)  S 
•  0 

O  J 

=  3,322 

L(N-1)+AL0G(FN) 
(-Jl+J2+J3)/2 
(+Jl-J2+J3)/2 
(+Jl+J2-J3)/2 
(-J1+L2+L3)/? 
(+Jl-L2+L3)/2 
(+Jl+L2-L3)/2 
(-Ll+J2+L3)/2 
(+Ll-J2+L3)/2 
(+L1+J2-L31/2 
=(-Ll+L2+J3)/2 
=(+Ll-L2+J3)/2 
=(+Ll+L2-J3)/2 
ED( 1)+MED(2)+MED(3) 
ED(4)+ME0(5)+MED(6) 
ED(7)+MED(8)+MED(9) 
rO( 10)+MED( 11)+ME0( 12) 
A( 1 )+MED( 12) 
A( 1)+MED( A) 
A( 1 )+MED(8) 

NE  MAXIMUM  OF  ( J  1+ J2+J 3) , ( J  1 +L2  +  L3 ) , < L 1  + J2  +  L3 ) , ( L 1+L2+J3 ) 
1) 

=  2,4 

-MAIN))  20,30,30 
N) 
E 

NE  MINIMUM  OF  ( J 1+ J2+L 1  +  L2 ) ,  ( J2+J3+L2+L3 ) , ( J3+J1+L3  +  L1 ) 
1) 

=  2,3 

-MB(N))  51,51,40 
) 
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MINP1=MIN+1 
MINI   =MINP1-MA(1) 
MIN2=MINP1-MA(2) 
MIN3=MINP1-MA(3) 
MIN4=MINP1-MA{4) 
MIN5=MINP1+1 
MIN6  =  Mti(  D-MIN 
MIN7=MB(2)-MIN 
MIN8=MB(3)-MIN 

SUM  SERIES  IN  OQURLE  PRECISION 
JK=10D-15 
S=l. 00-15 

IF  (KMAX)  65,65,55 
55  DO  60  K=1,KMAX 

UK=-UK*       DFLOAT(MINI-K)*DFL0AT(MIN2-K)*DFLOAT(MIN3-K)* 

1  DFL0AT<MIN4-K)/(DFL0AT(MIN5-K)*DFL0AT(MIN6+K.)*DFL0AT<MIN7+K)* 

2  DFLOAT (MIN8+K ) ) 

CUT  OFF  SERIES  AT  1.00-25 

IF(DABS(UK)-l.D-25)  65,60,60 
60  S=S+UK 
65  S=S*1.0D+15 

CALCULATE  DELTA  FUNCTIONS 

DELOG=0.0 

DO  70  N=l,12 

MUM=MED(M) 
70  DEL0G=DEL0G+FL(NUM+1) 

NUM1=MA( l)+2 

NUM2=MA(2)+2 

NUM3=MA(3)+2 

NUM4=MA(4)+2 

DEL0G=DEL0G-FL{NUM1)-FL(NUM2)-FL(NUP3)-FL(NUM4) 

DEL0G=G»5*DEL0G 

UL0G=FL(MIN5)-FL(MINI )-FL < MIN2 )-FL ( MIN3 ) -FL ( MIN4) -FL ( MIN6+1 )- 
1    FL(MIN7+1)-FL(MIN8+1) 

PL0G=D£L0G+UL0G 

IF(PL0G+64.)  72,75,75 

72  Q=PL0G+64. 
Q=EXP(Q) 
S6 J=w*S 

IF(ABS  (S6J)-1.)      73,73,74 

73  S6J=0. 
GO  TO  90 

74  S6J=S6J*EXP  (-64.) 
GO  TO  78 

75  P=EXP  (PLOG) 
S6J  =P*SNGL(S) 

78  MIN2=MIN/2 

IF  (MIN-2*MIN2)  80,90,80 


80 

S6J=-S6J 

90 

CONTINUE 

RETURN 

END 

c 

c 

SIBFTC  F3J 

FUNCTION 

LIST 

F3J{JD1,JD2,JD3,MD1,MD2,MD3) 
F3J  VERSION  I  FINAL  FORTRAN  IV 
DIMENSION  MTRU9) 
C0MM0N/FACT/FL(322) ,NCALL 
DOUBLE  PRECISION   UK,S,DFL0AT 
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50 
15 


1010 

1020 

1030 
1040 

1050 
1060 

1070 

30 
40 

45 

60 
70 
80 

90 

100 
110 
120 
130 


DFLOAT 

J1=JD1 

J2=JD2 

J3=J03 

M1=MD1 

M2=MD2 

M3=MD3 

IF(NCA 

NCALL= 

FL(1)= 

FL<2)= 

00    50 

FN=N-1 

FL(N)= 

1= Jl+J 

11=1/2 

IF{  1-2 

MTRH1 

I=J1-J 

11=1/2 

IF { 1-2 

MTRH2 

I=-J1+ 

11=1/2 

IF(I-2 

MTRK3 

IF(M1+ 

I=J1+M 

11=1/2 

IF(  1-2 

VITRU4 

MTRK5 

I=J2+M 

11=1/2 

IF ( 1-2 

MTRI (6 

MTRI (7 

1= J3+M 

11=1/2 

IF (  1-2 

MTRI (8 

MTRI (9 

DO    30 

IF(MTR 

CONTIM 

IF     ( J3 

KMIM=- 

GO    TO 

KMIN=0 

IF(-J3 

KMIN=- 

KMIN= 

IF( J2- 

KMAX=J 

30    TO 

KMAX=J 

IF(  J2  + 

KMAX=J 

KMAX=K 

MINI 


(  I)  =  I 


LL+1867)5,15,5 

-1867 
0.0 
0„0 
N=3,322 


FL(N-1 )+AL 
2-J3 

*I1)     1000, 
)  =  I1 
2+J3 

»I1)     1000, 
)  =  I1 
J2+J3 

*I1)     1000, 

)=I1 

M2+M3)     100 
1 

*U)     1000, 

)  =  I1 

)={ Jl-Ml)/ 

2 

*I1)     1000, 

)  =  I1 

)=( J2-M2)/ 

3 

*I1)     1000, 

)  =  I1 

)=( J3-M3)/ 

N=l,9 

I(N))     1000 

UE 

-J2+M1)     4^ 

J3+J2-M1 

60 


OG(FN) 

1010,1000 

1020,1000 

1030,1000 
0,1040,1000 

1050,1000 
2 

1060,1003 
2 

1070,1000 
2 

,30,30 
,45,45 


IN)  80,80,70 
100, IOC 


+J1+M2  -KM 

J3+J1+M2 

KMIN/2 

J3+M1)  90, 

1+J2-J3 

110 

1-M1 

M2-KMAX)  120,130,130 

2  +  M2 

MAX/ 2 

=  MTRI  (  D-KMIN+1 
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MIN2  =  MTRK5)-KMiN+l 

MIN3=MTRI (6)-KMIN+l 
MIN4=< J3-J2+Nl)/2+KMIN 
MIN5=( J3-Jl-N2)/2+KMIN 
UK  =loD-10 
S=1.D-10 
NCUT=0 

KNAX=KMAX-KNIN 
IF(KMAX)  165,  165,155 
155  DC  160  K=1,KMAX 

UK=-UK*DFLOAT(ylM       -K ) *C FLOAT ( MIN2-K } »DFLO AT ( M IN3-K ) /    {DFLOATtKMI 
1N+K)*OFLOAT(MN4  +  K)*CFLCAT{MIN5+K)) 
IF(DABS(UK)-1.D3G)  158,157,157 

157  UK=loD-10«UK 
S=1.0-10*S 
NCUT=NCUT+1 

158  IF (DABS(UK)-1„D-2C)  165,160,160 
16G  S=S+UK 

165  DEL0G=G.O 

DC     170    N=l,9 
NUM=MTRI (N) 

170  0EL0G=CELG6+FL(NUM+1) 
NUM={Jl+J2+J3)/2+2 
DEL0G=G35*(DELCG-FL(NUM)) 

ULCG=-FL(KMIK+1)-FL(MINI   )-F L ( MIN2 )-FL ( MIN3 )-FL< V IN4+1 )-FL ( NIN5+1 
1) 
PLCG=DELCG+ULCG 
IF(PLCG+80.0)  172,171,171 

171  IF(NCUT)  175,175,172 

172  SS=S 

SIC-  =SIGN(1,C,SS) 

SS=ABS(SS) 

SLCG=ALCG(SS)+FLCAT(NCUT+1)*AL0G( 1.E1C) 

F3J=SIG  *EXP(SLOG+PLCG) 

GC    T0178 
175    S  =  S*UC10 

P=EXP(PLOG ) 

F3J=P*SNGL(S) 
178  NUM=KMIM+( Jl-J2-F3)/2 

IF (M0C{NUf,2)  )  180,190,180 
180  F3J=-F3J 
190  CCNTINUE 

GC  TO  2000 
1C00  F3J=0.0 
2C00  RETURN 
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Listing  of  Version  II  Fortran  IV  Programs 

$I8FTC  F9J 

FUNCTION  F9J{ JDl , JD2 , JD3 , JD4 , JD 5, JD6, JD7 , JD8, JD9 ) 
C      F9J  VERSION  II  CALLS  S6J   FORTRAN  IV 

DIMENSION  MTRIA( 18)  ,KN(6)  ,KX(6) , NN ( 6) 

DOUBLE  PRECISION  SUM, S IG, FLK , S6 J , TERM ,DFLOAT 

DFLOAT( I )  =  I 

J1=J01 

J2=JQ2 

J3=JD3 

J4=J04 

J5=JD5 

J6=JD6 

J7=JD7 

J8=JD8 

J9=JD9 
C      ANGULAR  MOMENTUM  COUPLING  TESTS  FOR  9J  COEFFICIENT 

I=-J1+J2+J3 

11=1/2 

IF(  1-2*1  1 }  1000,1010,1000 
1010  MTRIA1 1 )  =  I1 

I=+J1-J2+J3 

11=1/2 

IF (  I-?* 1 1 )  1000,1020,1000 
1020  MTRIA(2)=I1 

I=+J1+J2-J3 

11=1/2 

IF (  1-2*1 1 )  1000,1030,1000 
1030  MTRIA(3)=I1 

I=-J4+J5+J6 

11=1/2 

I F( 1-2*1 1 )  100C, 1040, 1000 
1040  MTRIA(4)=I1 

I=+J4-J5+J6 

11=1/2 

IF ( 1-2*1 1 )  1000,1050,1000 
1050  MTRIA(5)=I1 

I=+J4+J5-J6 

11=1/2 

IF (  1-2*1 1 )  1000,1060,1000 
1060  MTRIA(6)=I1 

I=-J7+Jl3  +  J9 

11=1/2 

IF (1-2*11)  1000,1070,1000 
1070  MTRIA(7)=I1 

I=+J7-J8+J9 

11=1/2 

I  F (  1-2*11)  1000,1080,1000 
1080  MTRIA(b)=Il 

I=+J7+J3-J9 

11=1/2 

IF(I-2*U)  1000,1090,1000 
1090  MTRIA(9)=I1 

I=-J1+J4+J7 

11=1/2 

IF( 1-2* 1 1 )  1000,1100,1000 
1100  MTRIA(10)=I1 

I=+J1-J4+J7 

11=1/2 

IF (  1-2*1 1 )  1000,1110,1000 
1110  MTRIA(11)=I1 
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1120 


1130 


1140 


1150 


1160 


1170 


1180 


30 


35 


40 

50 


52 


I=+J1+ 
11=1/2 

IF(  1-2 

tfTRlA( 

I=-J2+ 

11=1/2 

IF(  1-2 

MTRlAt 

I=+J2- 

11=1/2 

IF(  1-2 

KTRIA( 

I=+J2+ 

11=1/2 

IF(I-2 

KTRIA< 

I=-J3+ 

11=1/2 

1F{  1-2 

M1RIA( 

I=+J3- 

11=1/2 

IF(  1-2 

MTRIA( 

I=+J3+ 

11=1/2 

IF(  1-2 

KTRIA( 

DO  30 

IF(MTR 

CONTl'-l 

KN(1)= 

KN(2 )= 

KN(3)= 

KN(4)= 

KN(5)= 

<N(6)= 

KX(1)= 

KX(2  )  = 

KX(3)= 

KX(4)= 

KX(5)= 

KX(6)= 

DO  35 

NN(K)= 

KSIGN= 

I=MIN0 

DO  40 

IFU-'M 

CONTIN 

KMIN=K 

KMAX=K 

GO  T0( 

J  =  Jl 

J1=J5 

J5  =  J 

J=J3 

J3=J8 

J8=J 

J=J6 


J4-J7 

*I1)  1000, 
12  )  =  1 1 
J5  +  J8 

*I1)  1000, 
1 3 )  =  1 1 

J5  +  J8 

*I1)  1000, 

14)=I1 

J5-J8 

*I1)  1000, 
15)=I1 
J6  +  J9 

*I1)  1000, 
16 ) =11 
J6  +  J9 


1120,1000 


1130,1000 


1140,1000 


1150,1000 


1160,1000 


*I1)  1000,1170,1003 

17 )  =  1  1 

J6-J9 


*I1) 
18)  =  I 
N=lt  1 
IA(N) 

UE 

MAXO( 

MAXOI 

MAXO< 

MAXO{ 

MAXO( 

MAXO( 

MINOl 

MINO( 

M  I  MO  ( 

MIMO( 

MIIMOI 

MINO  ( 

K=l,6 

KX(K) 

1 

(IMN(  1 

K=l,6 

N  (  K  )  ) 

U£ 

N  (  K  )  + 

X(K)  + 

130,5 


1000,  1183,1000 

1 

8 

)  i one, 30, 30 


I ABS( J2-J6) ,IABS(J1-J9),IABS(J4-J8)) 
IABS( J2-J7),  IA8S( J5-J9),  IABS( J4-J3) ) 
IABS( J6-J7) ,  IABS( J5-J1) , l\QSi J8-J3) ) 
IABS( J6-J1),  IABS( J2-J9),  IABS( J5-J7) ) 
IA8S( J2-J4),  IA8S(J3-J7),IA3S(J6-J8)) 
IABS( J3-J5),  IABS( J1-J8),  IA3S( J4-J9)  ) 
J2+J6, J1+J9, J4+J8) 
J2+J7, J5+J9, J4+J3) 
J6+J7, J5+J1, J8+J3) 
J1+J6, J2+J9, J5+J7) 
J2+J4, J3+J7, J6+J8) 
J3+J5, J1+J8, J4+J9) 

-KN(K ) 

)  ,  NN  (  2 )  ,  NN ( 3 ) , NN ( 4 ) , NN ( 5 ) , NN ( 6  )  ) 

40,50,40 

1 
1 
2,53, 54,55,56)  ,K 


52 


J6  =  J7 

J7  =  J 

GO  TO  130 

53 

J=J2 

J2  =  J7 

J7=J 

J=  J3 

J3  =  J4 

J4  =  J 

J=J5 

J5=J9 

J9  =  J 

GO  TU  130 

54 

J  =  J1 

J1=J2 

J2=J 

J=J4 

J4=J5 

J5  =  J 

J=J7 

J7  =  J8 

J8  =  J 

GO  TO  120 

55 

J=J1 

J1=J3 

J3  =  J 

J=  J4 

J4  =  Jo 

J6  =  J 

J=J7 

J7=J9 

J9  =  J 

GO  TO  12G 

56 

J  =  J2 

J2  =  J3 

J3=J 

J=J5 

J5  =  J6 

J6  =  J 

J  =  J8 

J8=J9 

J9  =  J 

120 

KSIGN=« 1-M0D( Jl+J2+J3+J4+J5+J6+J7+J8+J9,4) ) 

C 

SUMMATION  OF  SERIES  OF  EQUATION  (2)' 

130 

SUM=0„DO 

SIG=(-1 )**(KMIN-1)*KSIGN 

FLK=0FL0AT(KMIN) 

DO  200  K  =  KMIN,KN>AX,2 

TERM=FLK«S6J( Jl, J4»J7,J8, J9,K-1  )  *S6 J ( J2 , J 5, J 

1 

L*S6J( J3, J6, J9,K-l,Jl, J 2  ) 

FLK=FLK+2.D0 

200 

SUM=SUM+TERM 

F9J=SUM*SIG 

GO  TO  2000 

1000 

F9J=0«0 

2000 

RETURNS 

END 

C 

c 

$IBFTC  F6J      LIST 

,J4, K-ltJ6) 


53 


1000 
1010 

1020 

1030 

1040 

1050 

106C 

1070 

1080 

1090 

1100 

1110 


FUNCTI 

VERSIO 

DIMENS 

DOUBLE 

J1=JD1 

J2=JD2 

J3=JD3 

L1=L01 

L2=LD2 

L3=L03 

ANGULA 

I=-Jl+ 

11=1/2 

IF  (  I- 

F6J  =  0 

GO  TO 

MEO(l) 

I=+J1- 

11=1/2 

IF  (I- 

MED(2) 

I=+J1+ 

11=1/2 

IF  (  I- 

ME0(3) 

I=-J1+ 

11=1/2 

IF  (  I- 

MED(4) 

I=+J1- 

11=1/2 

IF  (  I- 

MED(d) 

I=+J1+ 

11=1/2 

IF  (  I- 

MED(6) 

I=-L1+ 

11=1/2 

IF  (  I- 

MED<7) 

I=+L1- 

11=1/2 

IF  (  I- 

MED(8) 

I=+L1+ 

11=1/2 

IF  (  I- 

MED(9) 

I=-L1+ 

11=1/2 

IF  (I- 

MEDI10 

I=+L1- 

11=1/2 

IF  (  I- 

MED( 11 

I=+L1+ 

11=1/2 

IF  (  I- 


ON  F6J{ J01,J02,J03,LD1?LD2,LD3) 
M  II  F6J  FUNCTION  CALLS  S6J   FORTRAN  IV 
ION  MEDU2) 
PRECISION   A6J     ,  S6J 


K  MOMENTUM  COUPLING  TESTS  FOR  6J  COEFFICIENT 
J2  +  J3 

2*11)  1000,1010,1000 

0 

100 

=  11 

J2  +  J3 

2*11)  1000,1020,1000 

=  11 

J2-J3 

2*11)  1000,1030,1000 
=  11 
L2  +  L3 

2*11)  100^,1040,1000 
=  11 
L2  +  L3 

2*11)  1000,1050,1000 

=  11 

L2-L3 

2*11)  1000,1060,1000 
=  11 
J2  +  L3 

2*11)  1000,1070,1000 
=  11 
J2  +  L3 

2*11)  1000,1080,1000 

=  11 

J2-L3 

2*11)  1000,1090,1000 
=  11 
L2  +  J3 

2*11)  1000,1100,1000 
)  =  I1 
L2  +  J3 

2*11)  1000,1110,1000 

)  =  I1 

L2-J3 

2*11)  1000,1120,1000 

54 


1120  MED(12)=U 

DO  ID  N=l,12 
IF  (MED(N))  1000,10,10 
10  CONTINUE 

F6J=S6J(J1,J2,J3,L1,L2,L3) 
100  RETURN 
END 


50 
15 


IBFTC  S6J 

DOUBLE 
VERSION 
COMMON/ 
DIMENSI 
DOUBLE 
1DSIGN,D 
DFLOAT( 
J1=JD1 
J2=JD2 
J3=J03 
L1=LD1 
L2=LD2 
L3=LD3 
DETERMI 
IF(NCAL 
5  NCALL=- 
CALCULA 
FL(1 )=0 
FL(2)=C 
DO  50  N 
FN=DFLO 
FL(N)=F 
MED(1)= 
MED<2)= 
MED(3)= 
MED(4) = 
MED(5)= 
MEO(b)= 
MED( 7)= 
MED(b)= 
MED (9)= 
MED(IO) 
MEDdl  ) 
MED(12  ) 
MA(  1  )=M 
MA(2)=M 
MA(3)=M 
MA(4)=M 
MB(1 )=M 
MB (2 )=M 
MB(3)=M 
DETERMI 
MAX=MA( 
DO  30  N 
IF  (MAX 
MAX=MA( 
CONTINU 
DETEKMI 
MIN=MB{ 
DO  51  N 


LIST 
PRECISIOfJ  FUNCTION  Sb J ( JD1 , J  02 , JD3» LD1 , LD2 , L03 ) 

II   FORTRAN  IV 
FACT/FLI322) , NCALL 
ON  MA(4),M3( 3) ,MEO{  12) 

PRECISION  FL,FN,A,OLOG,UK,S,OFLOAT,DELOG,ULOG,PLOG,SIG, 
ABS,SLOG,DEXP,P,Q 
I)=I 


20 
30 


E  WH 
+  186 
B67 
E  FL 
ODO 
ODO 

3,3 
T(N- 
(N-l 
-J1  + 
+J1- 
+  J1  + 
-J1  + 
+  J1- 
+  J1  + 
-L1  + 
+  L1- 
+  L1  + 
(-L1 
(+L1 
(+L1 
DID 
0(4) 
D(7) 
DUO 
(1)  + 
(1)  + 
(1)  + 
E  MA 
) 

2,4 
MA(N 
) 


ETHER  TO  CALCULATE  FL(N) 
7)  5,15,5 

(N)  S 


22 

1) 

J+DLOGl 

J2  +  J3), 

J2+J3) 

J2-J3) 

L2+L3), 

L2+L3) 

L2-L3) 

J2+L3) 

J2+L3) 

J2-L3). 

+L2+J3 

-L2+J3 

+L2-J3 

+MED(2 

+ME0I5. 

+ME0(8 

)+MED( 

MED( 12] 

MED( 4) 

MED( 8) 

XI  MUM  OF  ( J1  +  J2  +  J3) , ( J1  +  L2+L3),  (L1  +  J2  +  L3), (L1  +  L2+J3) 


))  20f30»30 


■  (FN) 
/2 
/2 
/2 
/2 
/2 
</2 
12 
/2 
/2 
)/2 
)/2 
)/2 

)+MED(3) 
)+MED(6) 
)+MED(9) 
U)+MED(  12) 
) 


E  MINIMUM 

) 

2,3 


OF  ( J1+J2+L1+L2) ,  ( J2+J3+L2+L3) , ( J3+J1+L3+L1) 


55 


IF  (MIN-MB(N))  51,51,40 
40  MIN=MB(N) 
51  CONTINUE 

KMAX=MIN-MAX 

MINP1=MIN+1 

MINI   =MINP1-MA(1) 

MIN2=MINP1-MA(2) 

MIN3=MINP1-MA<3) 

MIN4=MINP1-MA(4) 

MIN5=MINPL+1 

MIN6=MB(1)-MIN 

MIN7=MB(2)-MIN 

MIN8=MB(3)-MIN 
C      SUM  SERIES  IN  DOUBLE  PRECISION 

UK=1.D-15 

S=l. 00-15 

IF  (KMAX)  65,65,55 
55  DO  60  K=1,KMAX 

UK=-UK*DFL0AT{MINI-K)*DFLOAT(MIN2-K)*DFL0AT<MIN3-K)*DFL0AT(MIN4-K! 
l/CDFLOAT(MIN5-K)*DFLOAT(MIN6+K)*DFL0AT(MIN7+K)*DFL0AT(MIN8+K) ) 
C      CUT  OFF  SERIES  AT  1.0E-25 

IF(DABS(UK)-l.D-25)  65,65,60 
60  S=S+UK 
65  S=S»1. 00+15 
C      CALCULATE  DELTA  FUNCTIONS 

DEL0G=0.0D0 

DO  70  N=l,12 

NUM=MED(N) 
70  DEL0G=DELOG+FL(NUM+l) 

NUMl=MA(l)+2 

NUM2=MA(2)+2 

NUM3=MA(3)+2 

NUM4=MAI4)+2 

DEL0G=DEL0G-FL(NUM1)-FL(NUM2)-FL(NUM3)-FL<NUM4) 

DELOG=0.5D0«DELOG 

UL0G=FL(MIN5)-FL(MINI   )-FL( MIN2)-FL< MIN3)-FL (MIN4)-FL( MIN6+1)-  Ft 
1(MIN7+1)-FL(MIN8+1) 

PLOG=DELOG+ULOG 

IF(SNGL(PLOG)+64.0)  72,75,75 

72  Q=PL0G+64.D0 
Q=DEXP(Q) 
S6J=Q*S 
IF(DABS(S6J)-1.D0)   73,73,74 

73  S6J=0.D0 
SO  TO  90 

74  S6J=S6J*DEXPJ-64.D0) 
GO  TO  78 

75  P=DEXP(PLOG) 
S6J  =P*S 

78  MIN2=MIN/2 

IF  (MIN-2*MIN2)  80,90,80 
80  S6J=-S6J 
90  CONTINUE 
RETURN 
END 
C 
C 
ilBFTC  F3J      LIST 

FUNCTION  F3J(JD1,JD2,JD3,MD1,MD2,MD3) 
C      F3J  VERSION  II   FORTRAN  IV 

56 


C0MM0N/FACT/FH322)  ,NCALL 

DIMENSION  MTRK9) 

DOUBLE  PRECISION  FL,FN, A,DLOG»UK,S, DFLOAT, DELOG,ULOG ,PLOG,SIG, 
IDS IGN, DABS, S LOG, DEXP,P 

DFLOAT(I)=I 

J1=JD1 

J2=JD2 

J3=JD3 

M1=MD1 

M2=MD2 

M3=MD3 
;      DETERMINE  WHETHER  TO  CALCULATE  FL(N)  S 

IF(NCALL+1867)5,15,5 
5    NCALL=-1867 

FL(l)=O.DO 

FL(2)=0.D0 

DO    50    N=3,322 

FN=N-1 
50  FL(N)=FL(N-1)+DL0G(FN) 
15  I=J1+J2-J3 

11=1/2 

IF(I-2*I1)  1000,1010,1000 
1010  MTRI(1)=I1 

I=J1-J2+J3 

11=1/2 

IF(I-2*I1)  1000,1020,1000 
1020  MTRI(2)=I1 

I=-J1+J2+J3 

11=1/2 

IFC 1-2*11)  1000,1030,1000 
1030  MTRI(3)=I1 

IFIM1+M2+M3)  1000,1040,1000 
1040  I=J1+M1 

11=1/2 

IF ( 1-2*1 1)     1000,1050,1000 
1050    MTRK4)  =  I1 

MTRI(5)=(Jl-Ml)/2 

I=J2+M2 

11=1/2 

IF(I-2*I1)    1000,1060,1000 
1060    MTRI(6)=I1 

MTRI(7)=( J2-M2)/2 

I=J3+M3 

11=1/2 

IF(I-2*I1)    1000,1070,1000 
1070    MTRI(8)=I1 

MTRK9)  =  <J3-M3)/2 

DO    30    N=l,9 

IFIMTRKN})     1000,30,30 
30    CONTINUE 

IF    (J3-J2+M1)    40,45,45 
40    KMIN=-J3+J2-M1 

SO    TO    60 
45    KMIN=0 

60    IF(-J3+J1+M2    -KMIN)     80,80,70 
70    KMIN=-J3+J1+M2 
80    KMIN=KMIN/2 

IFIJ2-J3+M1)     90,100,100 
90    KMAX=J1+J2-J3 

GO    TO    110 
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100 
110 
120 
130 


155 


157 


158 
160 

165 


170 


171 
172 


175 


178 

180 

190 

1000 
2C00 


KM 
IF 
KM 
KM 
MI 
MI 
MI 
MI 
MI 
SU 
UK 
S= 
MC 
KM 
IF 
DO 
UK 
1N+ 
IF 
UK 
S  = 
NC 
IF 
S  = 
CA 
DE 
DO 
NU 
DE 
NU 
DE 
UL 
PL 
IF 
IF 
SI 
S  = 
SL 
F3 

go 

S  = 

P  = 
F3 
NU 
IF 
F3 
CO 
GO 
F3 
RE 
EN 


AX  =  J 
(  J2  + 
AX  =  J 
AX=K 
NI 

N2  =  M 
N3  =  M 
N4=( 
N5=( 
M  SE 

=  1. 
l.D- 
UT=0 
AX  =  K 
(KMA 

160 
=  -UK 
K)  *D 
(DAB 
=  1.0 
1«0— 
UT=^ 
(DAB 
S+UK 
LCUL 
LOG  = 

170 
M=MT 
LOG  = 
M=(J 
LOG= 
0G  =  - 
OG=D 
(SNG 
(NCU 
G  =0 
DABS 
OG  =  D 
J  =  SI 

TU1 
S*l. 
DEXP 
J-P* 
M=KM 
(MUD 
J=-F 
NT  1  J 

TO 

J  =  0e 

TURN 
D 


1-M 
M2- 
2  +  M 
MAX 
=  MT 
TRI 
TRI 
J3- 
J3- 
RIE 
D-l 
10 


,X)  120,130,130 


{  1  J-KMIN+1 

J-KMIN+1 

l-KMIN+1 

+M1J/2+KMIN 
-M2 1/2+KMIN 
IN  DOUBLE  PRECISION 


N 


MAX-KMI 
X)  165, 
K=l ,KM 
*DFLOAT 
FLOATIM 
S(UK)-1 
-10*UK 
1C»S 
CUT  +  1 
S(UK)-1.0-?' ) 165, 160, 16C 


165,155 
AX 

(MINI   -K)*DFL0AT(MIN2-K)*DFL0AT(MIN3-K )/ 
IN4+K)*DFL0AT(MIN5+K ) ) 
.030)158,157, 157 


(OFLC)AT(KMI 


ATE  DE 
Go  DO 
N=l,9 

run) 

DELOG+ 
1+J2+J 

oa5D0* 

FLIKMI 

ELOG+U 

L( PLOG 

T)  175 

SIGN! 1 

(S) 

LOG(S) 

G  *DEX 

78 

0+10 

(PLOG) 

S 

IN+U1 

(NUM,2 

3J 

ue 

2000 

u 


LTA  FUNCTIONS 


FL(NUM+1) 

3)/2+2 

(DELOG-FL(NUM) ) 

N+l )-FL(MINI )-FL( MIN2)-FL(M  IN3 )-FL ( MIN4+1 )-FL( MIN5+1 ) 

LOG 

)+80„0)     172,171,171 

,175,172 

.03, S) 


+DFLCAT(NCUT+1 ) *DLOG( 1. 0+1 0 ) 
P(SLOG+PLOG) 


-J2-M3J/2 

))     180,190,180 
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Listing   of    Fortran    IV   Testing    Programs 
SIBFTC    TST9J         LIST 
C  TEST    OF    F9J    FUNCTION    SUBPROGRAM       FORTRAN    IV 

DIMENSION  AJ(9) 
C      RUN  LABELLING, USE  FIRST  DATA  CARD  COLUMNS  1-72 
5  READ  (5,10) 
10  FORMAT(72H 

1  ) 

WRITE  <6,10) 
C      READ  MODE,  =1  DO  COEFFICIENTS  ONLY,  =  2  DO  BOTH  COEFFICIENTS  AND 
C      ORTHONORMALITY,  =  3  DO  ORTHONORMALITY  ONLY 
READ  (5,20)MODE 
IF(M0DE-2)  28,28,52 
C      READ  NMAX, NUMBER  OF  INDIVIDUAL  COEFFICIENTS  TO  BE  CALCULATED 
28  READ  (5,20)NMAX 

WRITE  (6,30)NMAX 
20  FORMAT(7I10) 

30  F0RMAT(24H  NUMBER  OF  CALCULATIONS=, 13  ) 
32  WRITE  (6,35) 

35  F0RMATI80H       Jl       J2       J3       J4       J5       J6       J7 
1       J8       J9       9J) 
DO  50  N=1,NMAX 
C      READ  TRUE  Jl , J2, J3, J4, J5, J6, J7, J8 , J9  FLOATING  POINT  E.G.  3/2=1.5 
READ  (5,40)(AJ(I ), 1=1,9) 
40  F0RMAT(9F8.1) 

Jl=2.0*AJ(l)+.001 
J2=2.0»AJ(2)+.001 
J3  =  2.0*AJ(3)+.Or>l 
J4=2.0*AJ(4)+.001 
J5=2.0*AJ(5)+.001 
J6=2.0»AJ(6)+.001 
J7=2.0*AJ(7)+.001 
J3=2.0*AJ(8)+.Or)l 
J9=2.O*AJ(9)+.O01 

ANSWER=F9J( J1,J2,J3,J4, J5,J6,J7,J8, J9) 
WRITE  (6,45) (AJ( I ), 1=1, 9), ANSWER 
45  F0RMAT(9F8.1,E15.7) 
50  CONTINUE 

52  IF(M0DE-2)    5,55,55 
C      9J  ORTHONORMALITY  CHECK 
C      READ  NMAX,  NUMBER  OF  NORMALIZATION  RUNS  TO  BE  MADE 

55  READ  (5,20)NMAX 
WRITE  <6,30)NMAX 
DO  100  N=1,NMAX 
WRITE  (6,56) 

56  FORMAT(72H       Jl       J2       J3       J4      J13    J13PR      J24 
1    J23PR        J  ) 

C      READ  TRUE  Jl , J2, J3, J4, J  13 , J 13PR , J24, J24PR, J    FLOATING  POINT 
READ  (5,40)TJ1,TJ2,TJ3,TJ4,TJ13,TJ13PR,TJ24,TJ24PR,TJ 
WRITE  (6,40)TJ1,TJ2,TJ3,TJ4,TJ13,TJ13PR,TJ24,TJ24PR,TJ 
J12MIN=2.0*ABS(TJ1-TJ2)+.001 
J12MAX=2.0*(TJ1+TJ2)+.001 
J34MIN=2.0*ABS(TJ3-TJ4)+.001 
J34MAX=2.0*(TJ3+TJ4)+.001 
J1=2.0*TJ1+.001 
J2=2.0*TJ2+.001 
J3=2»0*TJ3+.001 
J4=2.0*TJ4+.001 
J13=2.0*TJ13+.001 
J13PR=2.0»TJ13PR+.001 
J24=2.0»TJ24+.001 
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J24PR=2.0*TJ24PR+.031 
J=2.0*TJ+.001 
C      9J  ORTHONORMALITY  CHECK  BY  EQUATION  (4) 
SUM=0.0 

DO  100  J12=J12MIN,J12MAX,2 
DO  100  J34=J34MIN,J34MAX,2 
IF( J12+J34-J)  100,70,70 
70  IF(IABS( J12-J34)-J)80,80,100 

80  TERM=FLOAT< J12+1 ) *FLOAT ( J34+1 )*FLOAT( J13+1 )*FLOAT < J24+1 )  »  F9J(J1, 
1J2,J12,J3,J4,J34,J13,J24, J)*F9J(J1, J2, J  12, J  3, J4, J34, J13PR, J24PR, J) 
SUM=SUM+TERM 
TJ12=0.5*FLOAT( J12) 
TJ34=0.5*FL0AT(J34) 
WRITE  (6,60)TJ12,TJ34,TERM,SUM 
60  F0RMAT(2F6.1,2E15.7) 
100  CONTINUE 
GO  TO  5 
END 
C 

c 

SIBFTC  TST6J    LIST 

C      PROGRAM  FOR  TESTING  F6J  FUNCTION  SUBPROGRAM   FORTRAN  IV 

DIMENSION  AJ(6) 
C      RUN  LABELLING, USE  FIRST  DATA  CARD  COLUMNS  1-72 
5  READ  (5,10) 
10  F0RMAT(72H 

1  ) 

WRITE  (6,10) 
C      READ  MODE,  =1  DO  COEFFICIENTS  ONLY,  =2  DO  BOTH  COEFFICIENTS  AND 
C      ORTHONORMALITY,  =3  DO  ORTHONORMALITY  ONLY 
READ  (5,20)M0DE 
IF(M00E-2)  28,28,52 
C      READ  NMAX, NUMBER  OF  INDIVIDUAL  COEFFICIENTS  TO  8E  CALCULATED 
28  READ  (5,20)NMAX 

WRITE  (6,30)NMAX 
20  F0RMAT(7I10) 

30  F0RMAT(24H  NUMBER  OF  CALCUL AT IONS  =  ,  13  ) 
32  WRITE  (6,35) 

35  F0RMAT(55H       Jl       J2       J3       LI       L2       L3      6J) 
DO  50  N=1,,NMAX 
C      READ  TRUE  Jl , J2, J3, LI , L2,L3  IN  FLOATING  POINT  FORM  E.G.  3/2=1.5 
READ  (5,40)(AJ(I ), 1=1,6) 
40  F0RMAT(9F8.1) 

J1  =  2.0*AJ(  D+.001 

J2=2.0*AJ(2)+.001 

J3=2„0*AJ(3)+.001 

J4=2o0*AJ(4)+.001 

J5=2„0*AJ(5)+.001 

J6=2.0*AJ(6)+,001 

ANSWER=F6J(J1,J2,J3,J4,J5,J6) 

WRITE  (6,45) (AJ( I ), 1=1, 6), ANSWER 


ORTHOGONALITY  RUNS  TO  BE  MADE 


45 

F0RMAT(6F8.1,E15.8) 

50 

CONTINUE 

52 

IF(M0DE-2)    5,55,55 

c 

ORTHONORMALITY  CHECK 

c 

READ  IMAX,  NUMBER  OF 

55 

READ  (5,20)IMAX 
WRITE  (6,30)IMAX 
DO  150  1=1, IMAX 
WRITE  (6,56) 
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56  F0RMAT(48H       Jl       J2       J3       LI       L2     J3PR) 
C      READ  TRUE  J  1 ,  J2, J3, LI ,L2, J3PRIME  FLOATING  POINT  AS  ABOVE 
READ  (5,40)TJ1,TJ2,TJ3,TL1,TL2, TJ3PR 
WRITE  (6,40)TJ1,TJ2,TJ3,TL1,TL2,TJ3PR 
Jl=2.0*TJl+.001 
J2=2.0*TJ2+.001 
J3=2.0*TJ3+.001 
L1=2.0»TL1+.001 
L2  =2.0*TL2  +.001 
J3PR=2.C*TJ3PR+.001 

IF(IABS( J2-L1)-IABS( J1-L2) )  80,80,90 
80  NMIN=IABS( Jl-L2)+1 

GO    TO    100 
90    NMIN=IA3S(  J2-LD+1 
100    IF(J2         +L1         -Jl         -L2  )     110,110,120 

110    NMAX=J2  +L1  +1 

GO  TO  130 
120  NMAX=J1    +L2    +1 
C      6J  ORTHONORMALITY  BY  EQUATION  (3) 
130  SUM=0.0 

DO  150  N=NMIN,NMAX,2 

TERM=FLOAT(N)»FLOAT( J3PR  +1)*F6J(J1,J2,J3,L1,L2,N-1)»F6J(J1,J2,  J3 
1PR,L1,L2,N-1) 
SUM=SUM+TERM 
140  F0RMAT(3H  N=,I3,6H  TERM= , E16. 8, 5H  SUM=,E16.8) 

WRITE  <6,140)N,TERM,SUM 
150  CONTINUE 
GO  TO  5 
END 
C 
C 

$IBFTC  TST3J    LIST 

C      PROGRAM  FOR  TESTING  F3J  FUNCTION  SUBPROGRAM   FORTRAN  IV 
DIMENSION  AJ(3),AM(3) 
5  READ  (5,10) 
10  F0RMAT(72H 

1  > 

WRITE  (6,10) 
READ  (5,20)M0DE 
WRITE  (6,25)M0DE 
25  F0RMAT(6H  M0DE=,I1) 
IF(M0DE-2)  28,28,52 
28  READ  (5,20)NMAX 

WRITE  (6,30)NMAX 
20  F0RMAT(7I10) 

30  F0RMAT(24H  NUMBER  OF  CALCULAT IONS= , 13  ) 
32  WRITE  (6,35) 

35  F0RMAT(55H       Jl       J2       J3       Ml       M2       M3      3J) 
DO  50  N=1,NMAX 

READ  15,40) (AJJ I ),I=1,3),<AM( I), I =1,3) 
4C  F0RMAT(9F8.1) 

Jl=2.0»AJ(l)+.O01 
J2=2.0*AJ(2)+.001 
J3=2.0*AJ(3)+.001 
Ml=2.001*AM( 1) 
M2=2.001*AM(2) 
M3=2.001*AM(3) 

ANSWER=F3J( Jl , J2 , J3 , Ml , M2, M3 ) 

WRITE  (6,45) (AJ( I), I  =  1,3), (AM (I  ) , I  =  1 , 3) , ANSWER 
45  F0RMAT(6F8.1,E15.8) 
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50 
52 

55 


56 


140 
145 
150 


M3 


M3PR) 


CONTINUE 

IF(M0DE-2)    5,55,55 

ORTHONORMALITY  CHECK 

READ  (5,20)IMAX 

WRITE  (6,30)IMAX 

DO  150  I=1,IMAX 

WRITE  (6,56) 

F0RMATI48H       Jl       J2         J3    J3PR 

READ  (5,40)TJ1,TJ2,TJ3,TJ3PR,TM3,TM3PR 

WRITE  (6,40)TJ1,TJ2,TJ3,TJ3PR,TM3,TM3PR 

J1=2.0001*TJ1 

J2=2.0001*TJ2 

J3=2.0G01*TJ3 

J3PR=2.0001*TJ3PR 

M3=2»0001*TM3 

M3PR=2.0001*TM3PR 

SUM=0»0 

M1ST0P=2*J1+1 

M2STOP=2*J2+l 

DO  150  M=1,M1ST0P,2 

DO  150  N=1,M2ST0P,2 

M1=M-J1-1 

M2=N-J2-1 

TERM=FL0AT( J3+1)»F3J( Jl , J2, J3,M1 , M2 ,M3) *F3J ( Jl , J2 , J3PR , Ml ,M2 , M3PR) 

SUM=SUM+TERM 

IF(TERM-l.E-lO)  150,150,145 

F0RMAT14H  M1=I3,4H  M2=,I3,6H 

WRITE  <6,140)M1,M2,  TERM, SUM 

CONTINUE 

GO  TO  5 

END 


TERM=, E16.8,5H  SUM=,E16.8) 
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Appendix  III.   Sample  Input  for  Testing  Programs 
All  three  testing  (TST9J,  TST6J  and  TST3J)  programs  are  written 
in  a  similar  manner,  with  an  option  to  calculate  a  number  of  individual 
coefficients  only,  or  to  make  a  certain  number  of  orthonormality  runs 
only,  or  to  do  both. 


CARD 
NUMBER 


FORTRAN 
FORMAT 


EXPLANATION 


Next  NMAX 
cards 


Next  card 


Next  NMAX 
cards 


72H 


110 


110 


9F8.1 


110 


Run- labelling  card  used  for  entering  name, 
date,  specific  purpose  of  run,  remarks,  etc. 

MODE.   Mode  choice  for  test.   M0DE=1,  calculate 
individual  coefficients  listed  below  only; 
M0DE=2,  calculate  both  individual  coefficients 
and  perform  orthonormality  checks  as  well; 
M0DE=3,  do  orthonormality  checks  only. 

NMAX.      If   M0DE=1   or   =2,    NMAX=number   of   indi- 
vidual coefficients   to  be  calculated;   if 
M0DE=3,    NMAX=number   of   orthonormality  rela- 
tions  to  be   evaluated. 

If  M0DE=1  or  =2,  enter  arguments  of  the  indi- 
vidual coefficients;  if  M0DE=3,  enter  arguments 
of  the  orthonormality  relations.   (See  below 
for  argument  identification). 

NMAX.  This  card  appears  only  if  M0DE=2,  and 
gives  the  number  of  orthonormality  relations 
to  be  calculated. 


9F8.1      Arguments  of  the  orthonormality  relations. 

This  card  is  used  only  if  M0DE=2.   (See  below 
for  argument  identification. ) 
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Arguments  of  Individual  Coefficients: 

PROGRAM         Arguments  in  order  Notation 

TST9J  jl  »  JS  >  J3  >  J4  >  J5  >  JS  )  J7  »  J8  >  J9 


{ji    ja    J3^| 
J4    J5    J6  y 
J7     J8     J9j 


TST6J  J1.j2.j3  ,h>^^3  See  Eq.  (l')    in   text, 

TST3J  ji  ,  js  ,  J3  ,mi  ,ms  ,1113  See  Eq.  (6)    in  text. 
Arguments   for   Orthonormality   Relations: 

PROGRAM  Arguments    in   order  Notation 

TST9J  JI.JS  »J3  »J4,»J13»J13»J34,  ja.4  '»j  See  Ecl'  (4)     in  text ' 

TST6J  ja.js  »J3  *h>^s»33  '  See  Eq-  <3)  in  text. 

TST3J  ]i  ,  j3  ,  j3  ',1113  ,m3  See  Eq.  (7)  in  text. 
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5 

.5 

2.5 

2. 

3.5 

1.5 

2. 

4. 

3. 

1. 

.5 

3.0 

2. 

3.5 

1.5 

2. 

4. 

3. 

1. 

.5 

2.5 

2. 

3.5 

1. 

2. 

4. 

3. 

1. 

.5 

2.5 

2. 

3.5 

1.5 

2. 

5. 

3, 

2. 

1.5 

1.5 

0. 

3.5 

3.5 

4. 

4. 

2. 

4. 

40.       5.      RC.      42.      45.      45.      40.      39.      35. 
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2 

3 

8.0 

6.0 

5.0 

5.0 

6.0 

6.0 

8.0 

6.0 

5.0 

5.5 

4.5 

3.5 

8.0 

6.0 

5.0 

5.5 

4.5 

4.5 

2 

20. 

20. 

22. 

20. 

20. 

22. 

40. 

42. 

38. 

40. 

37. 

37. 
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2 

4 

8. 

8. 

8„ 

0. 

-8. 

8. 

8. 

8. 

8. 

-1. 

6. 

7. 

8. 

8. 

7. 

7. 

-8. 

1. 

8. 

8. 

6. 

8. 

-6. 

-2. 

2 

40. 

40. 

4C. 

40. 

0. 

0. 

80. 

80. 

80. 

80. 

0. 

0. 
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